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ABSTRACT 


'J  The  Electrical  and  Computer  Engineering  Department  at  the  Naval  Postgrad¬ 
uate  School  has  a  need  for  additional  software  to  be  used  in  instructing  students 

•4.1“-, 

studyingisignal  processing.  This  software  will  be  used  in  a  PC  lab  or  at  home. 

This  thesis  provides  a  set  of  disks  written  in  APL  (A  Programming  Language) 
which  allows  the  user  to  input  arbitrary  signals  from  a  disk,  to  perform  various 
signal  processing  operations,  to  plot  the  results,  and  to  save  them  without  the  need 
for  complicated  programming. 

The  software  is  in  the  form  of  a  digital  signal  processing  ^toolkit.”  The  user 
can  select  functions  which  can  operate  on  the  signals  and  interactively  apply  them 

in  any  order.  The  user  can  also  easily  develop  new  functions  and  include  them  in 

?•  <> 

the  “toolkit.” 

The  thesis  includes  brief  discussions  about  the  library  workspaces,  a  user  man¬ 
ual,  function  listings  with  examples  of  their  use,  and  an  application  paper.  The 
software  is  modular  and  can  be  expanded  by  adding  additional  sets  of  functions. 
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I.  INTRODUCTION 


The  Electrical  and  Computer  Engineering  Department  at  the  Naval  Postgrad¬ 
uate  School  has  a  need  for  signal  processing  software  packages  for  use  by  students 
studying  digital  signal  processing  (DSP).  The  purpose  of  this  thesis  is  to  provide  a 
set  of  tools  in  the  form  of  APL  workspaces  on  floppy  disks  that  students  can  use  at 
home  or  at  school  on  IBM-PC-compatible  computers  for  solving  digital  signal  pro¬ 
cessing  problems.  Workspace  is  a  collection  of  functions  and  data  that  are  placed 
together  in  order  to  do  any  job. 

The  software  supports  the  following  currently  offered  courses  in  DSP: 

1.  EC  2400-Discrete  Systems 

2.  EC  3400-Introduction  to  Digital  Signal  Processing 

3.  EC  3410-Introduction  to  Discrete  Time  Random  Processes 

4.  EC  4430- Advanced  Signal  Processing  and  Spectral  Estimation 

5.  EC  4440-Multidimensional  Digital  Signal  Processing 

Specific  APL  workspaces  were  developed  for  EC3400,  EC3410,  and  EC4440. 
The  software  uses  STSC  Corporation’s  APL*P.uUS  PC  and  STATGRAPHICS  sys¬ 
tems  which  are  both  site  licensed  to  the  Naval  Postgraduate  School.  APL*PLUS 
PC  is  STSC’s  version  of  the  APL  language  for  the  IBM-PC  and  compatibles.  Sim¬ 
ilar  versions  of  the  APL  language  axe  also  available  from  STSC  for  other  com¬ 
puters.  STATGRAPHICS  (Statistical  Graphics  System)  is  a  PC  software  package 
integrating  a  variety  of  statistical  functions  with  color  graphics.  STATGRAPH¬ 
ICS  is  written  in  APL  (specifically  APL*PLUS  PC)  and  offers  direct  access  to  the 


host  language.  Currently  STATGRAPHICS  is  available  only  for  the  IBM-PC  and 
compatibles. 

The  DSP  software  developed  in  this  thesis  is  in  the  form  of  a  digital  signal 
processing  “toolkit.”  The  user  can  select  functions  to  operate  on  the  signals  and 
interactively  apply  them  in  any  order.  The  user  can  also  easily  develop  new  func¬ 
tions  and  include  them  in  the  “toolkit.”  The  tools  axe  provided  on  a  set  of  floppy 
disks  that  users  can  carry  to  their  own  PC  machine.  The  sets  of  disks  contain: 

1.  The  APL*PLUS  Executive  Program 

2.  STATGRAPHICS  disks 

3.  UTILITY  workspace 

4.  EC3400  workspace 

5.  EC3410  workspace 

6.  EC4440  workspace 

The  main  products  of  the  research  are  the  disks  that  include  the  signal  pro¬ 
cessing  library  workspace,  and  the  documentation  provided  in  this  thesis.  A  typical 
library  contains  functions  such  as:  Fast  Fourier  Transform  (FFT),  Inverse  FFT, 
Sine,  Cosine,  Sample  Autocorrelation,  Mean,  and  other  similar  functions.  In  de¬ 
veloping  the  DSP  software  the  following  requirements  were  met: 

1.  Ability  to  input/output  from/to  a  disk  data  file  in  free  format.  This  allows 
convenient  communication  with  DOS  and  with  programs  written  in  other  lan¬ 
guages  or  on  other  computer  systems. 

2.  Ability  to  provide  for  graphic  output  to  represent  signals,  frequency  responses, 
and  other  functions  including  2-D,  3-D,  and  contour  plotting. 

3.  Ability  to  perform  filtering,  convolution,  and  other  common  signal  processing 
operations  in  the  signal  (time/space)  domain. 


v.v.v. .  .  . 


4.  Ability  to  compute  transfer  functions,  to  compute  frequency  responses,  and 
generally  to  perform  complex  arithmetic  and  other  common  frequency  domain 
operations. 

5.  Ability  to  generate  random  signals. 

6.  Ability  to  compute  statistical  characteristics  of  the  signals  including  correla¬ 
tion  functions,  spectra,  and  prediction  errors. 

The  user  can  combine  functions  in  any  desired  order  to  solve  DSP  computer 
assignments  and  can  concentrate  on  the  signal  processing  exercises  without  the 
need  to  do  complicated  programming. 

Although  the  intent  was  to  develop  the  signal  processing  toolkit  specifically  for 
the  IBM-PC  and  compatibles,  the  system  is  portable.  STSC  has  APL  available  on 
IBM  mainframes,  Digital  Equipment  Corporation  VAX,  (under  UNIX  and  VMS), 
Apple  Macintosh,  and  others.  These  DSP  workspaces  axe  portable  to  all  of  these 
systems  with  few,  if  any,  changes.  In  addition,  the  DSP  workspaces  are  easily 
carried  over  to  other  APL  systems  with  minor  changes  to  only  some  of  the  most 
basic  functions;  bits  and  pieces  of  the  software  have  existed  for  some  time  under 
VS- APL  on  the  IBM  mainframe  and  under  UNIX  4.3  bsd. 

The  remainder  of  the  thesis  is  organized  as  follows. 

Chapter  II  contains  a  brief  description  of  the  APL  language  and  the  STAT- 
GRAPHICS  software  system.  Chapter  III  deals  with  the  DSP  library  and  discusses 
the  application  of  the  four  workspaces.  The  Utility  workspace  includes  general 
functions  such  as  GETDATA  and  PUTDATA  for  the  data  input/output  to  DOS; 
the  EC3400  workspace  includes  basic  one-dimensional  DSP  functions;  the  EC3410 
workspace  includes  one-dimensional  DSP  functions  for  statistical  signal  processing; 
and  the  EC4440  workspace  includes  multidimensional  DSP  functions.  Chapter  IV 
outlines  the  conclusions,  and  the  benefits  of  using  the  entire  package  as  a  software 
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toolkit.  Appendix  A  provides  a  simple  user  manual  for  the  software  package.  Ap¬ 
pendix  B  provides  a  function  library  with  a  brief  description  of  each  function.  Ap¬ 
pendix  C  contains  function  listings  and  samples  of  their  use.  Appendix  D  furnishes 
samples  of  computer  assignments  and  their  solutions  using  the  tools  developed  in 
this  thesis. 


II.  APL  AND  STATGRAPHICS  ENVIRONMENT 


The  signal  processing  software  is  written  in  APL  and  uses  the  APL*PLUS  PC 
interpreter  and  the  STATGRAPHICS  software  system.  Although  STATGRAPH¬ 
ICS  provides  many  statistical  functions  in  itself,  only  the  plotting  functions  from 
STATGRAPHICS  axe  used.  Following  is  a  brief  description  of  APL  and  the  STAT¬ 
GRAPHICS  system. 

A.  APL  [1] 

Originally  APL  began  as  a  notation  to  express  mathematical  procedures.  Its 
originator,  Dr.  Kenneth  Iverson,  published  his  notation  in  1962  as  A  Program¬ 
ming  Language  from  which  the  name  “APL”  is  derived  [2]. 

One  of  the  strongest  advantages  of  APL  is  its  ability  to  handle  complicated 
array  data-structures  without  extensive  programming.  The  language  allows  the 
definition  of  numerical  u  ms  or  a  series  of  items  with  one  symbol.  Examples  of 
data  that  can  be  represented  are: 

•  SCALAR  -  A  single  value 

•  VECTOR  -  A  group  of  scalars 

•  MATRIX  -  A  group  of  vectors 

•  ARRAY  -  A  group  of  matrices 

The  language  makes  it  as  easy  to  perform  an  arithmetic  operation  using  a 
matrix  or  vector  argument  as  to  perform  the  operation  using  simple  numbers. 
This  characteristic  is  illustrated  in  the  following  example. 


1 


M«-l6 

M 


(define  vectors  with  6  numbers) 


N«-2xi6 

N 

2  4  €  3  10  12 


n*m  (vector  multiplication) 

2  3  18  32  50  72  v 


M<-3  3  p<-9 
M 

12  3 
4  5  6 
7  8  9 


(define  matrices  with  3  rows  and  3  columns) 


N«-3  3  p3*l9 
N 

3  6  9 

12  15  18 
21  24  27 


NxM 

3  12  27 

48  75  108  (matrix  arithmetic  operation) 

147  192  243 


Notice  that  when  two  items  are  of  the  same  dimensions,  the  arithmetic  operation  is 
performed  element  by  element.  Thus  N  x  M  is  not  the  usual  linear  algebra  product 
of  matrices.  The  latter  is  realized  by  a  composite  operation  denoted  by  +.  x. 


Scalar  arguments  are  automatically  extended  to  be  conformable  with  the  arrays. 
Thus  the  statement  Mx2  yields 

2  4  6 

8  10  12 
14  16  18 

There  is  frequently  no  need  for  loops  or  similar  mechanisms  as  in  other  lan¬ 
guages  such  as  Fortran.  From  the  user’s  point  of  view  the  language  performs  the 
operation  on  the  elements  of  the  array  in  parallel  without  the  need  for  loops. 

Note  also  that  there  is  no  need  for  dimensions  because  APL  carries  the  dimen¬ 
sions  of  an  item  along  with  its  definitions. 

APL  provides  a  function  to  help  the  user  to  define  or  to  determine  the  number 
of  elements  in  the  variable.  The  function  represented  by  the  Greek  character  up” , 
is  called  “shape,”  since  it  gives  the  form  of  the  variable.  On  the  previous  page 
there  is  tin  example  where  we  define  the  matrices  M  and  N  using  the  function  “p” . 
In  order  to  find  the  shape  of  any  variable,  use  the  function  “p”  with  the  variable 
name  following  it.  The  following  is  an  example. 


APL  has  a  large  variety  of  mathematical  functions  that  can  easily  manipulate 
the  data  and  can  handle  operations  on  a  simple  as  well  as  on  a  complicated  level. 
Thus,  it  is  possible  to  place  on  a  single  line  a  combination  of  several  functions  which 
work  together  in  a  pipeline.  All  of  this  can  be  done  interactively.  The  flexibility  in 
manipulating  array  data  easily  makes  the  language  well  suited  for  single  processing 
analysis.  The  following  are  a  few  examples: 


(outer  product)  used  to  define  separable 
2-P  function  x(nl,n2)  =  xl(nl)x2(n2) 


2345  ».*  6  "1 
3  1  ~l  S 

9  2  0  7 

10  3  1  3 

11  4  2  9 


"3  4 


H«-3 

It 


3  pi  35790334 


1  3  5 
7  9  0 
3  5  4 


The  domino  operator  (calculate  M  inverse) 

3M 

"0.1313131313  "0.09090909091  0.2272727273 

0.1414141414  0. 1313131313  "0. 1767676763 

0.1515151515  "0.09090909091  0.06060606061 


cm  3  6 


(solve  Mx  —  C  for  x) 

"0.9090909091  1.04040404  1.757575758 

(check  solution  by- 
multiplying  constant  vector 

CSM3  xC  r  w  •  s 

"0.9090909091  1.04040404  1.757575753  ^  M  averse  j 

+.  x  is  matrix  product 


APL  uses  symbols  to  represent  its  own  built-in  functions  and,  in  addition,  APL 
allows  the  user  to  define  other  functions.  The  solution  to  one  problem  consists  of 
a  collection  of  the  programmer’s  functions,  each  one  of  which  can  be  applied  inter¬ 
actively  by  the  user  or  can  be  organized  to  call  another  function.  The  advantage 
of  this  organization  is  the  ease  of  debugging  and  that  only  a  few  lines  of  APL  can 
do  volumes  of  work. 

APL  has  a  concept  called  the  workspace  which  is  a  collection  of  functions  and 
data  that  are  placed  together  in  order  to  do  a  job.  That  collection  is  stored  in  the 
computer  memory  during  an  APL  session  and  gives  the  operator  immediate 
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feedback  and  accessibility,  and  the  workspace  can  be  saved  on  disk  and  loaded  back 
into  the  memory  at  any  time  for  future  use. 

B.  STATGRAPHICS  SOFTWARE  PACKAGE  [3] 

STATGRAPHICS  is  a  Statistical  Graphics  System  written  in  APL  for  scientific 
and  engineering  applications  and  for  data  analysis.  The  software  supports  dot 
matrix  printers,  monochrome/color  terminals,  and  pen  plotters.  It  runs  under  the 
APL*PLUS  PC  interpreter  version  5.0  or  later,  which  allows  the  user  to  write 
his  own  functions.  Using  STATGRAPHICS  in  this  mode  requires  familiarity  with 
APL. 

The  system  also  has  a  menu  that  contains  graphic  functions,  a  set  of  data 
analysis  and  statistical  procedures.  This  menu  is  illustrated  in  Figure  2.1.  In 
addition,  the  menu  contains  data  management  and  system  utilities  that  allow  the 
user  to  select  the  system  environment,  control  color,  plot  dimensions,  and  so  on. 
The  most  often  used  screens  for  DSP  applications  are: 

1.  Plotting  functions  (Figure  2.2) 

2.  Time  Series  Analysis  (Figure  2.3) 

3.  Experimental  Design  (Figure  2.4) 

STATGRAPHICS  has  options  for  modifying  graphs  and  for  on-line  help. 

To  use  STATGRAPHICS  with  APL  the  user  needs  to  start  APL*PLUS  PC  as 
usual  and  then  load  the  STATGRAPHICS  workspace  from  the  start-up  disk.  By 
going  back  and  forth  from  the  APL  environment  to  STATGRAPHICS,  the  user 
can  process  data  and  plot  the  results  easily.  This  makes  the  combination  of  APL 
and  STATGRAPHICS  powerful  and  very  convenient  to  use. 


I  - 

DATA  MANAGEMENT  AND  SYSTEM  UTILITIES 
A.  Data  Management 
a.  System  Environment 
c.  Sapors  writer  and  Graphics  Replay 

O.  Plotter  Interface 

I 

PLOTTING  AND  DESCRIPTIVE  STATISTICS 
E.  Plotting  Functions 

P.  Descriptive  Methods 

G.  Estimation  and  Testing 

H.  Distribution  Functions 

I.  Exploratory  Data  Analysis 


TIME  SERIES  PROCEDURES 

L.  Forecasting 

M.  Quality  Control 
M.  Smoothing 

O.  Time  Series  Analysis 

ADVANCED  PROCEDURES 

P.  Categorical  Data  Analysis 

Q.  Multivariate  Methods 

R.  Nonparametric  Methods 

S.  Sampling 

T.  Experimental  Design 


ANOVA  AND  REGRESSION  ANALYSIS 

J.  Analysis  at  Variance 

K.  Regression  Analysis 


MATHEMATICAL  AND  USER  PROCEDURES 

U.  Mathematical  Functions 

V.  Supplementary  Operations 


Use  cursor  Iceys  to  highlight  desired  section.  Then  press  ENTER. 
lHelp  2 Edit  3Savscr  4Prtscr  S  SGo  7Vars  8Cad  9Reviev  lOQuit 

INPUT  11/18/87  13:08  STATGRAPHICS  Vers.  2.1 


Figure  2.1.  STATGRAPHICS  Main  Menu. 


PLOTTING  FUNCTIONS 


1.  X-Y  Line  and  Scattarplots 

2.  Multiple  X-Y  Plots 

3.  X-Y-Z  Line  and  Scattarplots 

4.  Multiple  X-Y-Z  Plots 

5.  Barcharts 
8.  Piecharts 

7.  Component  Line  Charts 


Use  cursor  Keys  to  highlight  desired  procedure.  Then  press  ENTER. 

IHelp  2Edit  3Savscr  4Prtscr  5  SGo  7Vars  8Cad  9 Review  lOQuit 

INPUT  11/18/37  13 : 10  STATGRAPHICS  Vers.  2.1 


Figure  2.2.  The  Plotting  Function  Menu 


Ml 


TIME  SERIES  ANALYSIS 


Horizontal  Tina  Sequence  Plot 
Vertical  Tina  Sequence  Plot 
Saasonal  Subsarias  Plot 
Autocorrelation  Function 
Partial  Autocorralation  Function 
Cross -Correlation  Function 
Sinpla  or  Saasonal  Differencing 
Mean  or  Trend  Removal 
Box-Cox  Transformation 
Pariodogran 

Integrated  pariodogran 
Data  Tapering 

Plotting  vs.  Fourier  Frequencies 
Sox-JenJcins  ARIMA  Modeling 
Cross-Correlation  Matrix  Plot 


Use  cursor  keys  to  highlight  desired  procedure.  Then  press  ENTER. 
lHelp  2Edit  SSavscr  4Prtser  5  6 So  Wars  scad  9Review  lOQuit 

INPUT  11/18/87  13:13  STATGRAPHICS  Vers.  2.1 


Figure  2.3.  The  Time  Series  Analysis  Menu. 


EXPERIMENTAL  DESIGN 


1.  Full  and  Fractional  Factorials 

2.  Central  Composite  Designs 

3 .  Alias  structure 

4.  Response  Surface  Plotting 


Usa  cursor  keys  to  highlight  dasired  procedure.  Then  press  ENTER. 
lHelo  2Edit  3Savscr  4Prtscr  5  SGo  Wars  8Cmd  9Revi*v  toou’t 


Figure  2.4.  The  Experimental  Design  Menu. 


in.  THE  DSP  WORKSPACES  AND  APPLICATIONS 

The  following  is  the  principle  part  of  the  thesis  research.  This  chapter  describes 
the  workspaces  that  were  developed  by  the  author,  and  the  way  they  can  be  applied. 

A.  UTILITY  WORKSPACE 

This  workspace  contains  functions  of  general  use  that  can  be  applied  in  combi¬ 
nation  with  functions  from  other  workspaces  in  this  package.  Appendix  B  contains 
the  list  of  workspace  functions  and  a  brief  discussion  of  each  one. 

Two  of  the  most  useful  functions  in  the  utilities  are  described  and  demon¬ 
strated  below.  The  functions  GETDATA  and  PUTDATA  allow  the  user  to  transfer 
data  from  the  APL  environment  to  the  DOS  environment  and  vice  versa.  These 
functions  use  file  system  commands  that  are  specific  to  STSC’s  implementation  on 
the  IBM-  PC.  As  such,  the  exact  APL  code  for  GETDATA  and  PUTDATA  is  not 
always  portable  to  APL  on  other  systems.  However,  versions  of  GETDATA  and 
PUTDATA  can  be  written  for  any  system.  Versions  are  currently  available  for  the 
IBM  mainframe  under  VS/APL  and  for  a  version  of  APL  under  UNIX  4.3  bsd  on 
the  VAX.  To  the  user,  GETDATA  and  PUTDATA  provide  a  common  interface  for 
communication  with  other  programs. 

The  following  is  an  example  showing  the  transfer  of  data  using  the  GET¬ 
DATA  and  PUTDATA  functions.  Figure  3.1  describes  two  sets  of  data  in  an  APL 
workspace.  In  order  to  transfer  those  sets  to  a  file  in  the  DOS  environment  the 
user  needs  to  type  the  following  commands: 
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X  PUTDATA  *  X . DAT ’ 

Y  PUTDATA  ‘ Y . DAT  * 

X  is  5  x  5  and  Kisa2x4x4  array 


V, 

a 


£ 

31 


I 


. 

r 

* 

i 


0.923 

X 

1.35 

2.773 

3.7 

4.625 

5.53 

6.475 

7.4 

3.325 

9.25 

10.175 

11.1 

12. 025 

12. 93 

13.375 

14.3 

15.725 

16.65 

17.573 

13.5 

19.425 

20.35 

21.273 

22.2 

23.125 

X 

1260 

2260 

3260 

4260 

5E60 

6260 

7260 

3260 

9260 

1261 

1244 

1244 

"3244 

"4244 

"3244 

"6244 

"7244 

"3244 

"9244 

"1245 

"1.1245 

"1.2245 

“1.3245 

1.4245 

"1.3243 

"1.6245 

1.7245 

1.3245 

"1.9245 

"2245 

"2.1245 

1.2245 

Figure  3.1.  X  and  Y  Data  Sets. 


The  left  argument  to  PUTDATA  is  the  variable  name;  the  right  argument  is  the 
file  name  or  path.  Note  that  after  transferring  any  data  set  to  DOS  the  data  is  not 


I 


3 

tT. 

I 


i 

V 

V 


1 


erased  and  remains  available  in  the  workspace. 

When  terminating  APL,  the  user  is  placed  in  the  DOS  environment.  Figure 
3.2  shows  that  the  data  set  transfer  was  successfully  completed.  Note  that  the  first 
line  of  the  data  set  gives  the  dimension  or  “shape”  of  the  data. 

By  going  bank  to  APL  environment  and  typing  the  commands: 

X<— GETDATA  ‘X.DAT’ 

Y*— GETDATA  ‘Y.DAT’ 

the  data  is  transferred  correctly  from  the  DOS  environment  to  the  APL  workspace. 
Figure  3.3  shows  that  the  data  has  the  same  original  shape.  The  data  axe  read 
into  the  workspace  and  assigned  to  the  variables  X  and  Y .  The  data  still  remain 
available  as  a  file  in  the  DOS  environment. 
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l-iviv.v: , :  £ £ 


)OFT 


(terminate  APL) 


C:\APL>TYPE  X.DAT 

5  5 

0.925  1.35  2.775  3.7 

4.625  5.55  6.475  7.4 

3.325  9.25  10.175  11.1 

12.025  12.95  13.375  14.3 

15.725  16.65  17.575  13.5 

19.425  20.35  21.275  22.2 

23 . 125 

C:\APL>TYPE  Y . DAT 


2  4  4 


1E60 

2ES0 

3E60 

4E60 

5E60 

6E60 

7E60 

3E60 

9E60 

1E61 

-1E44 

-2E44 

-3E44 

-4E44 

-5E44 

-6E44 

-7E44 

-3E44 

-9E44 

-1E45 

-1.1E45 

-1.2E45 

-1.3E45 

-1.4E45 

-1.5E45 

-1.6E45 

-1.7E45 

-1.3E45 

-1.9E45 

-2E45  -2 

. 1E45  -2 

-2E45 

C:\APL> 


Figure  3.2.  X  And  Y  Data  Sets  In  DOS  Environment. 


<>x 

5  5 

X 

0.925  1. 35  2.775  3.7  4.625 

5.55  6.475  7.4  3.325  9.25 

10.175  11.1  12.025  12.95  13.375 

14.3  15.725  16.65  17.575  13.5 

19.425  20.35  21.275  22.2  23.125 


(type  of  contents  of  files) 


pY 

2  4  4 


Y 


1E60 

2ES0  3E60 

4E60 

5E60 

6E60  7E60 

3ES0 

9E60 

1E61  "1E44 

'2E44 

~3E44 

“4E44  ”5E44 

"5E44 

~7E44 

“3E44  “9E44 

"1E45 

-1.1E4S 

"1.2E45  "1.3E45 

“1.4E45 

"1.5E45 

"1.6E45  "1.7E45 

“1. 3E45 

“1.9E45 

”2E45  '2.1E45 

"2.2E45 

Figure  3.3.  Getting  X  And  Y  Into  APL  Workspace. 

B.  EC3400  WORKSPACE 

EC  3400  is  a  course  at  the  Naval  Postgraduate  School  that  introduces  first 


principles  of  digital  signed  processing.  The  topics  covered  include:  the  Discrete 


Fourier  Transform  and  the  FFT  algorithm,  flow-graph  and  matrix  representation 


of  filters,  ideal  filters  and  approximation,  and  design  of  recursive  and  nonrecursive 


digital  filters,  fast  convolution  and  correlation. 


The  following  workspace  enables  the  student  to  solve  computer  assignments  on 


the  above  subjects  in  a  simple  way.  This  workspace  contains  functions  such  as:  lin¬ 


ear  convolution,  circular  convolution,  FFT,  inverse  FFT,  and  complex  arithmetic. 


Some  of  those  functions  will  be  described  later.  In  addition,  the  EC3400  workspace 


contains  functions  which  allows  the  student  to  design  one- dimensional  filters.  The 


methods  used  are  the  Fourier  method  and  the  frequency  sampling  method  and 


windowing  of  the  filters  coefficients.  The  full  description  of  the  workspace  is  given 


in  Appendix  B.  The  following  is  an  example  of  how  to  use  this  workspace. 


EC  3400  Computer  Assignment  f4l 


A  highpass  filter  is  to  have  an  analog  cutoff  frequency  of  6kHz.  The  filter  is 


to  be  implemented  by  a  digital  filter  having  a  sampling  frequency  of  40kHz.  Plot 
the  frequency  response  of  the  filter  with  and  without  a  Hamming  window  if  51 


coefficients  are  used. 


Solution 


The  solution  given  here  uses  functions  provided  in  the  EC3400  workspace.  In 
an  actual  classroom  assignment,  students  may  be  asked  to  generate  their  own  func¬ 
tions  or  to  solve  the  problem  in  their  own  way.  By  using  the  function  DIGFREQ, 
the  digital  cutoff  frequency  is  calculated.  (Figure  3.4).  This  function,  like  several 
others  is  included  for  convenience  only;  the  calculation  in  converting  radian  to  dig¬ 
ital  frequency  can  be  done  easily  in  the  workspace  by  using  the  built-in  arithmetic 
operations.  Figure  3.5  describes  the  ideal  frequency  response  to  be  generated. 


vc 


By  using  the  function  HPCOEFF,  51  causal  impulse  response  coefficients  axe  cal¬ 
culated.  (See  Figures  3.6  and  3.7).  The  function  HPCOEFF  uses  the  following 
equation: 


h Lp{n )  =  k/ir(n  -  I)  sin([n  -  I]QC),  n  =  0, 1, 2, . . . ,  21 


TC«-40000  DIGFREQ  3000 
TC 

1.236637061 

TC+PI 

0.4 


Figure  3.4.  Calculating  a  Digital  Cutoff  Frequency. 


Figure  3.5.  The  Ideal  Frequency  Response. 
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HHJ^Sl  HPCOEET  TC 
hSP 

6. 3633 1673 3E "17  "0.01251377881  "3. 134689424E"3  8. 504448034E"3  0.01441574 
2. 33986119 1E"1 7  "0.01593319429  "0.01039432538  0.01100575628 
0.01392066822  2. 339861191E“17  "0.02162362082  "0.01439214283 
0.01559148806  0.02752097195  2.339861191£"17  "0.0336367435  "0.02338 
0.02672826525  0.05045511524  2.339861191E“17  "0.07563267236 
"0.06236595225  0.09354892838  0.3027306915  "0.6  0.3027306915 
0.09354392838  "0.06236595225  "0.07568267286  2. 339861191E"17 
0.05045511524  0.02672826525  "0.02338723209  "0.0336367435  2.3398611 
0.02752097195  0.01559148806  “0.01439214283  "0.02162362082 
2. 339861i91£"17  0.01392066822  0.01100575628  "0.01039432538 
*0.01593319429  2. 339861191E"17  0.01441574721  3 . 504448034E"3 
"8. 134689424E"3  "0.01261377831  6. 363316738E~17 

Figure  3.6.  Calculating  the  Impulse  Response  Coefficients. 


51  coeff.  for  causal  nonrecursive 

HP  filter 


Figure  3.7.  The  Impulse  Response  Coefficients. 


The  function  HAMMING  generates  “I71  samples  of  a  shifted  Hamming  window 


using  the  following  equation: 


W(n)  =  0.54  +  0.46  cos 


,0  <  n  <  I-  1 


This  is  multiplied  by  the  filter  coefficients  to  obtain  the  windowed  coefficients.  (See 


Figures  3.8  and  3.9). 

By  using  the  function  FREQRES,  the  frequency  response  with  and  without  a 
window  is  generated.  This  is  done  by  calculating  the  DFT  of  the  filter  coefficients 


i 
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W«- HAMMING  51 


W 

0.08087246878  0.08783237415  0.1016466798  0.122105975  0.1489001176  0.1816 
0.2197783849  0.2627880673  0.31  0.3606984983  0.4141150246  0.4694398 
0.5258342731  0.5824434454  0.6384092133  0.6928832078  0.7450396437 
0.7940878876  0.3392844181  0.3799441019  0.9154505797  0.9452656094 
0.9689372255  0.9861065906  0.9965134344  1  0.9965134344  0.9861065906 
0.9689372255  0.9452656094  0.9154505797  0.8799441019  0.8392844181 
0.7940878876  0.7450396437  0.6928832078  0.6384092183  0.5824434454 
0.5258342731  0.4694398388  0.4141150246  0.3606984983  0.31  0.2627880 
0.2197783849  0.1816229357  0.1489001176  0.122105975  0.1016466798 
0.08783237415  0.08087246878 


hHPW«-hHpxW 


Figure  3.8.  Calculating  the  Hamming  Window  Samples. 


O  lO  20  30  40  50  60 

samples  C  nl 


Figure  3.9.  The  Hamming  Window. 

and  finding  the  magnitude  of  the  result  as  shown  in  a  section  of  the  FREQRES 
function  (Figure  3.10).  The  following  describes  the  APL  commands  which  generate 
the  frequency  responses  and  Figures  3.11  and  3.12  show  the  corresponding  results. 

HHP<— 128  FREQRES  hHP 
HHPW<— 128  FREQRES  hHPW 
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VFREQRESCU]V 
HW«-M  FREQRES  h:OIO;N:H 
a 

A  GENERATING  FREQUENCY  RESPONSE  MAGNITUDE  AND  PHASE  FROM  TIME  DOMAIN 
A  COEFFICIENTS 

A 

010«-0 

N«-<M-ph)pO 

-■h-h.N  _ , 

H«-FFT  h 

HW«-CM+2)fCXMAGN  H) 

PW<-CM+2mXPHAS  H) 

•♦0 

A  THE  FUNCTIONS  FFT . XMAGN . XPHASE  ARE  USED. 

A  IN  ORDER  TO  RECEIVE  A  BETTER  RESOLUTION  THE  COEFFICIENTS  ARE  ZERO- 
A  PADDET  TO  'M'  SAMPLES. 

A  M*-NUMBER  OF  SAMPLES  INCLUDING  ZERO  PADDING  (USE  ONLY  RADIX  TWO  NO 
A  h«-THE  FILTER  COEFFICIENTS  C VECTOR] 

A  HW: FREQUENCY  RESPONSE  MAGNITUDE  CM+2  SAMPLES) 

A  PW: FREQUENCY  RESPONSE  PHASE  CM+2  SAMPLES)  CRAD. ] 


Figure  3.10.  Using  the  Function  FFT. 


£0]  CxPI/643 

Figure  3.11.  The  Frequency  Response  Using  51  Coefficients. 

C.  EC3410  WORKSPACE 

ECTMIO  is  a  course  that  gives  an  introduction  to  discrete-time  random  process. 
The  topics  covered  are:  description  of  discrete-time  random  signals  and  random 
vectors,  linear  transformations,  sampling  of  continuous- time  random  signals,  esti¬ 
mation,  and  spectral  analysis.  The  EC3410  workspace  was  developed  in  order  to 
enable  the  student  to  solve  computer  assignments  in  digital  signal  processing  and 


CxPI/643 


Figure  3.12.  The  Frequency  Response  Using  Hamming  Window. 

learn  to  write  algorithms  in  an  efficient  way  without  doing  excessive  programming. 
This  workspace  contains  functions  that  compute  the  mean,  circular  convolution, 
linear  convolution,  fast  Fourier  transform,  power  spectrum  and  other  quantities. 
A  full  description  of  the  workspace  functions  is  given,  in  Appendix  B. 

By  using  the  STATGRAPHICS  software  package  with  the  EC3410  workspace, 
a  graphic  capability  is  achieved.  The  following  is  a  computer  assignment  in  statis¬ 
tical  dsp  and  its  solution  using  the  EC3410  workspace  and  the  STATGRAPHICS 
software. 

Note  that  in  the  actual  assignment,  students  are  required  to  write  some  of  the 
functions  that  are  used  below.  However  developing  the  functions  in  APL  using  the 
tools  provided  is  not  difficult. 


EC  3410  COMPUTER  ASSIGNMENT  [5] 

1.  Three  data  sets  SLl.DAT,  SL2.DAT,  and  SL3.DAT  are  given  on  a  disk.  The 
mean  of  the  signal  data  SLl.DAT  is  0.75.  Subtract  this  constant  value  from 
the  signal  to  generate  a  new  zero-mean  signal  SLl'  that  you  will  use  in  this 
assignment.  The  other  data  sets  (SL2,  SL3)  have  zero  mean.  Do  the  remaining 
parts  for  each  of  the  data  sets  SLl',  SL2,  and  SL3. 

2.  Using  the  correlation  method,  generate  a  3  by  3  Toeplitz  correlation  matrix 
for  the  signal  data.  Print  this  matrix. 

3.  Solve  a  set  of  Normal  equations  involving  the  correlation  matrix  that  was  just 
generated  to  obtain  the  2nd  order  linear  predictive  filter  parameters  and  the 
prediction  error  variance. 

4.  Apply  the  filter  to  the  original  data  set  and  generate  the  prediction  error 
signal.  Plot  this  signal  and  compute  its  variance.  Does  it  compare  well  with 
the  theoretical  prediction  error  variance  you  obtained  by  solving  the  normal 
equations? 


Solution 

1.  By  using  the  function  MEAN,  the  SL1  data  set  mem  is  computed  and  is 
subtracted  from  the  original  SLl.  The  new  data  set  SL1'  has  a  zero-mean. 
(Figure  3.13) 
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Figure  3.13.  Using  the  Function  MEAN. 
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Figure  3.14.  Using  the  Functions  SACF  to  Evaluate  the  Covariance  for 
SLl.DAT. 
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Figure  3.15.  Using  the  Functions  SACF  and  COV  for  SL2  and  SL3. 

2.  By  using  the  function  SACF  (sample  autocorrelation  function),  the  first  three 
correlation  function  lags  axe  calculated  and  then  by  using  the  function  COV  a 
Toeplitz  matrix  is  generated  from  the  correlation  function  data.  The  process 
described  in  the  above  paragraph  is  repeated  on  the  SL2  and  SL3  data  sets, 
as  seen  in  Figures  3.14  and  3.15. 

3.  The  function  LPFP  (linear  prediction  filter  parameters)  solves  the  Normal 
equations  and  it  calculates  the  2nd  order  filter  coefficients  and  the  prediction 
error  as  demonstrated  in  Figure  3.16.  (This  could  also  be  done  without  the 
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Figure  3.16.  Solving  the  Normal  Equations  Using  the  Function  LPFP. 

LPFP  function  by  setting  up  simple  linear  equations  and  solving  them  as 
described  in  Chapter  II.) 

4.  Using  the  LCV  (linear  convolution)  function,  the  three  data  sets  axe  convolved 
with  the  2nd  order  filter  parameters  that  were  found  in  the  previous  step  to 
provide  the  prediction  error  sequence.  The  prediction  errors  of  the  signals 
are  determined  by  computing  the  initial  value  of  the  autocorrelation  function 
of  the  prediction  error  sequence  as  seen  in  Figure  3.17.  As  can  be  seen,  the 
results  are  very  similar  to  the  prediction  error  values  derived  by  the  analysis 
of  Figure  3.16.  Figure  3.18  is  a  plot  of  the  prediction  sequence  for  each  data 
set. 

Additional  examples  of  using  the  EC3410  workspace  axe  provided  in  Appendix 
D. 

D.  EC4440  WORKSPACE 

EC  4440  is  a  course  at  the  Naval  Postgraduate  School  in  multidimensional 
digital  signal  processing.  The  course  deals  with  the  analysis  of  signals  that  are 
functions  of  the  two  or  more  independent  variables.  The  course  develops  both 
time/space  and  frequency  domains  approaches  and  highlights  the  subjects  that 
are  different  from  one- dimensional  signal  processing.  Some  of  the  topics  covered 
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Figure  3.17.  Calculating  The  Prediction  Error  Using  The  Function  PRER. 

are:  two-dimensional  circular  and  linear  convolution,  difference  equations,  recur¬ 
sively  computable  systems,  Fourier  analysis,  and  methods  of  2-D  filter  design.  The 
EC4440  workspace  enables  the  student  to  solve  computer  assignments  in  multidi¬ 
mensional  signal  processing  by  using  functions  such  as:  2-D  convolution,  2-D  Fast 
Fourier  Transform  (FFT)  and  inverse  2-D  FFT,  and  filter  design  using  the  win¬ 
dowing  method  or  the  McClellan  transformation.  The  plots  are  provided  by  using, 
as  before,  the  STATGRAPHICS  software  package.  More  detailed  lists  and  descrip¬ 
tions  of  the  workspace  functions  axe  provided  in  Appendices  B  (DSP  library)  and 
C  (functions  listings  and  samples  of  runs). 
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The  following  is  an  example  of  how  to  design  a  lowpass  filter  using  the  workspace 
functions.  The  design  here  is  done  by  the  McClellan  transformation  method. 
Again,  as  in  other  examples,  the  students  are  asked  to  write  some  of  these  functions. 


COMPUTER  ASSIGNMENT  ON  McCLELLAN  TRANSFORMATION  [6] 
Consider  the  prototype  filter: 


—7 r  —  Air  Air  ir  w 

1.  Using  a  32  point  FFT,  find  the  impulse  response  hp(n)  and  the  coefficients 
a(n)  in  the  representation 


Hp(u>)  =  ^  a(n)Tn  [cos  w]  N  =  15 


2.  Now  let  F(u>  1,0/2)  =  \  (—1  +  cosu/i  +  cos u?2  +  cos u\  •  cos^) 

Compute  H(ui,u}2)  =  J^n=o  a(n)Tn[.F(u>i,u/2)]  and  generate  both  3-D  and 
contour  plots  of  this  frequency  response. 

Note:  You  can  take  advantage  of  the  recursion 

Tn(x)  =  2xTn-i(x)  -  r„_2(x) 

and  the  built-in  recursive  nature  of  APL  functions  to  generate  the  Chebyshev 
polynomials. 


H 

*3 


Solution 


The  function  PROFILT  is  used  to  specify  a  32-sample  ideal  lowpass  prototype 
filter;  then  by  using  the  function  COEFF  the  impulse  response  coefficients  are 
calculated.  As  shown  in  the  following  steps  the  COEFF  function  uses  the  IFFT 
function  to  calculate  the  impulse  response,  h(n),  of  the  prototype  filter. 


fcoeffco:* 

CO]  a«-COEFF  Hp:aiO:SI2E:HH:H 
Cl]  A 

C21  A  COMPUTING  THE  COEFFICIENTS  FROM  THE  PROTOTYPE  FILTER  C USING  IFFT 3 
C3 1  A 

C4  ]  aia«-o 

C33  SIZE«-(pHp3+2 
C6  7  a«-SriEaO 
C  7  ]  H-IFFT  Hp 


The  filter  coefficients  are  defined  as: 


,  \  A 

a(n)  = 


h(Q ),  n  =  0 
2 h(n),  n  >  0 


and  are  found  by  performing  the  following  steps  in  the  function  COEFF. 
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Figure  3.19.  Calculating  the  Impulse  Response  Coefficients. 
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Figure  3.20.  1-D  Prototype  Filter. 
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The  above  calculations  axe  demonstrated  in  Figure  3.19.  (Notice  that  the  prototype 
filter  response  is  defined  on  the  interval  0  to  2ir  rather  than  —  ir  to  +7r.)  Figures 
3.20  and  3.21  are  plots  of  the  prototype  filter  and  the  prototype  filter  coefficients. 

The  ideal  filter  would  have  a  step  discontinuance  at  the  points  u  =  OAir 
and  1.6tt  but  this  appears  in  the  plot  as  slanted  lines  because  the  plot  uses  finite 
numbers  of  samples. 
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Figure  3.21.  The  Prototype  Filter  Coefficients. 

The  function  TRANSFNC  generates  the  transformation  function  by  imple¬ 
menting  the  following  equation: 

F(w  1,  w2)  =  0.5(— 1  -j-  cos  wl  +  cos  w2  -f  cos  u;l  •  cos  iv2). 

The  surface  plot  capability  of  STATGRAPHICS  is  used  to  form  a  3-D  plot  and  a 
contour  plot  of  the  transformation  function  as  seen  in  Figure  3.22. 

By  using  the  function  MCCLEL  the  frequency  response  of  the  desired  filter  is 
generated.  The  function  generates  the  nth  Chebyshev  polynomial,  Tn[F(wl,  w2)) 
by  using  the  recursion  T„(x)  =  2 (x)T„_i(x)  —  Tn- 2(x).  By  using  the  following 


summation: 


^a(n)rn[F(a>l,u>2)] 


m7TnrrrrFrrTn7\ 


the  frequercy  response  H(wl,w2)  is  achieved.  The  above  steps  are  shown  in  a 
section  of  the  MCCLEL  function. 
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Figure  3.23  shows  the  frequency  response  plot  of  the  lowpass  filter.  In  order 
to  apply  a  status  report  during  calculation,  a  status  printout  was  added  to  the 
function  as  seen  in  Figure  3.24. 


E.  DISCUSSION 

The  workspaces  and  functions  developed  in  this  thesis  provide  a  complete 
set  of  tools  for  solving  the  problems  currently  assigned  to  students  in  courses  EC 
3400,  EC  3410,  and  EC  4400.  They  also  provide  capabilities  to  solve  many  other 
problems  in  DSP. 
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Figure  3.24.  Status  Printout  During  Calculation  of  Chebyshev  Polyno¬ 
mial. 

Some  of  the  functions  developed  and  included  in  the  workspace  are  finite 
specialized  and  were  specifically  coded  as  examples  of  solving  the  particular  DSP 
homework  assignments.  Others  are  provided  as  tools  that  can  be  applied  to  many 
different  problems.  As  already  mentioned,  it  is  very  easy  to  expand  this  set  of  tools 
to  include  others. 

In  using  the  workspaces  and  functions,  the  instructor  will  usually  provide 
only  a  subset  of  the  functions  developed  here  to  students  enrolled  in  the  course. 
Students  would  then  use  the  functions  provided,  as  well  as  the  general  features  of 
APL  to  solve  the  assigned  problems.  Since  the  workspaces  are  extremely  modular, 
an  instructor  will  find  it  easy  to  select  only  those  functions  that  he  wants  students 
to  have  and  provide  others  at  a  later  time. 


IV.  SUMMARY  AND  CONCLUSION 


The  purpose  of  this  thesis  was  to  develop  a  convenient  set  of  software  tools  for 
students  and  instructors  in  order  to  solve  computer  assignments  in  digital  signal 
processing  on  personal  computers  using  APL. 

The  main  requirement  was  achieved  by  developing  APL  workspaces  with  sets  of 
functions  that  could  be  applied  to  perform  DSP  operations.  Computer  assignments 
in  EC  3400.  EC  3410,  and  EC  4440  are  easily  solved  by  using  this  software  package. 
The  user  needs  to  have  basic  knowledge  in  APL  and  needs  to  know  how  to  use  part 
of  the  STATGRAPHICS  package  in  the  APL  environment.  However,  this  basic 
knowledge  is  easily  acquired. 

The  structure  of  the  APL  software  is  very  flexible  and  convenient  to  use.  By 
combining  a  few  functions,  the  user  can  achieve  complicated  calculations  with  little 
programming  effort.  Any  workspace  can  be  expanded  by  the  student  or  professor 
by  adding  new  functions. 

The  results  of  this  thesis  demonstrate  that  the  APL  is  easy  and  convenient 
to  use,  especially  in  solving  digital  signal  processing  problems.  This  thesis  will 
hopefully  make  it  easier  for  new  students  to  start  using  the  language. 


APPENDIX  A 
USER  MANUAL 


A.  SYSTEM  REQUIREMENTS 

The  software  package  contains  the  following  disks: 

1.  APL*PLUS  Executive  Program  version  5.1  or  later  by  STSC 

2.  STATGRAPHICS  -  Statistical  Graphics  System,  version  2.1  or  later  by  STSC 

3.  Utility  workspace,  and  EC3410,  EC4440,  and  EC3400  workspaces 

STATGRAPHICS  requires  the  following  hardware  and  system  software: 

1.  IBM  PC,  PC-XT,  PC- AT  or  compatible 

2.  512K  RAM 

3.  Keyboard  with  APL  characters 

4.  Two  double-sided,  double-density  disk  drives  or  one  floppy  disk  and  hard  disk. 

5.  Graphics  adapter  card 

6.  DOS,  version  2.0  or  later 

7.  Dot  matrix  graphics  printer  (optional) 

S.  Math  copi  ocessor  (optional  but  recommended) 

The  software  package  can  be  operated  from  floppy  disks  or  from  the  hard  disk. 

B.  PREREQUISITES  FOR  USER 

In  order  to  use  the  package  efficiently,  the  user  should  have  a  basic  knowledge 
of  APL. 

In  addition,  it  is  strongly  recommended  that  a  copy  of  the  STATGRAPHICS 
manual  [3]  be  available  to  the  user  since  this  Appendix  refers  the  user  to  Ref.  3 
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for  more  explanations.  The  user  should  not  have  loaded  any  special  DOS  resident 
programs  since  this  may  interfere  with  the  proper  execution  of  STATGRAPHICS 
and  cause  the  system  to  “hang  up.” 

C.  HOW  TO  GET  INTO  THE  DSP  LIBRARY 

Follow  the  instructions  below  to  start  the  session: 

1.  Start  the  computer  using  DOS. 

2.  If  you  are  using  a  dual  floppy  disk  drive,  place  the  APL*PLUS  executive  disk 
in  the  default  disk  drive.  If  you  axe  using  a  hard  disk  system,  make  sure  that 
the  directory  with  the  APL  and  STATGRAPHICS  program  (usually  the  APL 
directory)  is  your  current  directory. 

Enter:  APLCOM  if  you  are  using  APL  version  5.1  or 
Enter:  A  if  you  are  using  APL  version  6.4. 

These  commands  call  up  batch  files  specifically  for  each  version.  A  batch  file 
contains  a  set  of  commands  that  must  be  executed  before  it  runs  the  APL  program 
in  order  to  have  the  software  matched  with  the  specific  hardware. 

The  system  will  respond  as  follows: 


I  APUPLU3  PC  APPLICATION  DEVELOPMENT  SYSTEM  I 

I - I 

I  Version  5.1  Serial  Number  134257  Copyright  1985,  ST3C,  Inc.  I 

I  All  rights  reserved.  Unauthorized  reproduction  of  this  software  I 

I  is  prohibited  and  violates  U.S.  Copyright  Laws.  APL*PLU3  is  a  I 

I  registered  trademark  and  service  mark  of  STSC,  Inc.  I 

*  .  —  .  ■  ■  -  —  ... — - - 

CLEAH  US  11/20/1987  14:31:54 

3.  Place  the  desired  DSP  library  disk  in  the  default  disk  drive,  as  for  example, 
the  disk  containing  the  UTILITY  workspace. 

Enter:  )LOAD  UTILITY 

If  you  are  using  a  hard  disk  and  want  to  load  the  workspace  from  drive  A, 
enter:  )LOAD  0  UTILITY.  The  “0”  (zero)  stands  for  drive  A,  “1”  stands  for 
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drive  B,  and  “2”  stands  for  drive  C  (hard  disk).  The  system  will  respond  as 
follows: 

JLOAD  0  UTILITY 

0  UTILITY  SAVED  12/13/1987  13:12:28 


NAVAL  POSTGRADUATE  SCHOOL 


MONTEREY . CALIFORNIA 


VERSION  1.0 


UTILITIES  WORKSPACE 

THESIS  BY 
Y.KATZIR.I.A.F. 

ADVISOR: PROF.  C.».  THERRIEN. 

SEPTEMBER  1987 


At  this  point  you  can  use  the  UTILITY  workspace  as  described  in  Appendices 
B,  C,  and  D. 

4.  Before  terminating  the  session  save  your  results  by  entering: 

)SAVE  UTILITY  or  )SAVE  0  WSID  (WSID  means  your  file  name) 
the  system  will  respond: 

)SAVE 

0  UTILITY  SAVED  11/19/1987  00:  24:06 

5.  To  terminate  APL  enter:  )OFF  and  you  will  be  automatically  returned  to 
DOS. 


D.  USING  STATGRAPHICS  WITH  APL 
Follow  the  instructions  below  to  start  the  session: 

1.  Start  the  computer  under  DOS.  If  a  graphics  printer  is  attached  be  sure  that 
the  printer  is  set  up  for  graphics  printing.  (DOS  has  to  run  graphics.com  when 
system  boots.) 

2.  Load  APL  as  described  in  the  previous  section  (C.b.) 

3.  After  the  APL*PLUS  PC  system  is  loaded,  replace  the  APL  disk  with  the 
STATGRAPHICS  start-up  disk  and  enter: 


)LOAD  STATGRAF 


ii 

C; 


WARNING:  When  using  the  STATGRAF  workspace  use  only  the  )PCOPY 
command  to  enter  functions  or  variables  into  the  workspace. 


During  the  loading  procedure,  the  system  will  display  a  copyright  banner  and 
the  following  message: 

“SYSTEM  INITIALIZING.  PLEASE  BE  PATIENT.  THIS  WILL  TAKE  A  FEW  MOMENTS. 

When  initialization  is  complete,  the  system  will  prompt  for  questions  that  axe 
related  to  the  computer  configuration.  Select  the  appropriate  answers  (See  Ref. 

3,  pg.  3-1  to  3-3  for  more  details.)  The  system  will  ask  if  you  want  to  save  the 
settings  previously  entered  in  order  to  have  an  automatic  logon  the  next  time  you 
start  STATGRAPHICS.  The  Main  Menu  will  then  be  displayed  as  in  Figure  A.l. 
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J.  Analysis  of  variance 

K.  Regression  Analysis 


TIME  SERIES  PROCEDURES 
L.  Forecasting 
N.  Quality  Control 

N.  Smoothing 

O.  Tims  Series  Analysis 

ADVANCED  PROCEDURES 

P.  Categorical  Data  Analysis 

Q.  Multivariate  Methods 

R.  Nonparametric  Methods 

S.  Sampling 

T.  Experimental  Design 

MATHEMATICAL  AND  USER  PROCEDURES 

U.  Mathematical  Functions 

V.  Supplementary  operations 


Use  cursor  leeys  to  highlight  desired  section.  Then  pree#  ENTER. 
lHslp  2 Edit  SSevscr  4Prtscr  5  SGo  7Vars  3 Cad  9Review  lOQuit 

INPUT  11/18/87  13:08  STATGRAPHICS  Vers.  2.1 


Figure  A.l.  Statgraphics  Main  Menu. 


4.  When  you  axe  in  the  Main  Menu,  press  the  ESC  key  to  enter  APL  immediate 

execution  mode.  Execute  the  function.  Enter:  RESTART  to  return  to  the 

menu-control  mode. 

Since  the  complete  set  of  functions  in  the  STATGRAPHICS  system  exceeds 
the  memory  space  available  on  most  personal  computers,  you  should  load  only 
those  functions  and  variables  that  are  necessary  to  perform  the  desired  tasks.  For 
DSP  applications  the  major  procedures  used  are  described  in  the  following  section. 

D.  STATGRAPHICS  MAJOR  PROCEDURES  USED 

After  starting  STATGRAPHICS,  the  first  display  is  the  Main  Menu  (See  Fig¬ 
ure  4..1.) 

The  Main  Menu  contains  22  sections  under  six  major  categories. 

The  Main  Menu  gives  you  access  to  all  the  data  analysis  and  graphics  proce¬ 
dures. 

A  highlighted  cursor  bar  is  displayed  on  the  Main  Menu  screen.  In  order  to 
select  any  desired  Menu  Procedure,  move  the  cursor  using  the  arrow  keys  and  press: 
ENTER.  STATGRAPHICS  procedures  axe  stored  in  four  disks.  When  using  floppy 
disks  carefully  follow  the  instruction  on  the  screen. 

The  following  is  a  short  description  of  the  most  frequently  used  procedures  in 
this  DSP  software  package: 

1.  System  Environments  [3] 

The  STATGRAPHICS  system  environment  section  is  under  category: 

“DATA  MANAGEMENT  AND  SYSTEM  UTILITIES” 

This  section  allows  you  to  modify  the  STATGRAPHICS  system  environ¬ 
ment  to  fit  your  needs.  This  section  has  nine  procedures.  Those  most  often  used 
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Blinking  CO-1): 


Figure  A. 2.  The  Screen  Options  Menu. 

1.  Screen  option:  Allows  you  to  change  colors,  text,  and  properties  of  the  graphics 
screen.  (See  Figure  A.2). 

2.  Graphic  options:  Allows  you  to  change  sizes,  orientation,  and  positions  of  the 
graphics  screen.  (See  Figure  A. 3.) 
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Figure  A.3.  Graphics  Options  Menu. 


2.  Plotting  Functions  [3] 

Plotting  functions  are  found  on  the  Main  Menu  under  the  category: 
“PLOTTING  AND  DESCRIPTIVE  STATISTICS” 
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STATGRAPHICS  provides  seven  plotting  procedures  you  can  use  to  pro¬ 
duce  a  variety  of  plots.  The  following  is  an  example  of  using  the  X-Y  Time  Plot 


procedure. 


1.  Start  in  the  APL  environment  and  define  the  variables  X  and  Y  as  follows: 

0  STATGRAP  APL  Ins 

You  are  now  under  control  of  the  APL  interpreter. 

To  resume  menu  control,  type: 

RESTART 

- *  X«-t.lO 

X 

123456789  IQ 

- *  Y«-l  2  2.5  3.3  4  5  6  7.5  8.3  9.5 

Y 

1  2  2.5  3.3  4  5  6  7.5  3.3  9.5 


Enter:  RESTART 


Select  the  X-Y  Line  and  Scatterplots  procedures  from  the  Main  Menu. 


2.  Fill  in  the  Menu  entry  as  demonstrated  in  Figure  A. 4. 


(USE  CAPITAL  LETTERS) 


'  X:  X 
►  Y:  Y 


X-Y  Line  and  Scatterplots 


Log 

Mo 

Mo 


Point  labels: 

Point  codes: 
Point  colors: 
Std.  err.  X: 
Std.  err.  Y: 


Points:  Yes 
’  Lines:  Yes 


Figure  A.4.  X-Y  Plotting  Menu 

3.  Press  the  F6  key  in  order  to  see  the  plot.  (See  Figure  A. 5.) 

4.  Press  the  F4  key  to  print  the  screen,  (if  graphics  printer  is  attached) 

5.  Press  the  F5  key  to  make  any  changes  in  the  plot  (including  titles). 

6.  Press  the  F10  key  to  leave  the  screen  and  to  return  to  the  menu. 


.•»  h.ii.ii.i  ■■■  u».i ’.Mit.m.M.m.v.*! 
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Figure  A.5.  The  Plot  of  Y  vs  X. 

3.  Response  Surface  Plotting  [3] 

This  procedure  produces  3-D  surface  plots  and  contour  plots,  for  previ¬ 
ously  generated  two  dimensional  functions.  The  procedure  is  under  section:  “Ex¬ 
perimental  Design”  on  the  main  menu  under 

“ADVANCED  PROCEDURES” 

The  following  is  an  example  of  how  to  use  this  procedure: 

1.  Start  with  a  definition  of  a  matrix  of  values  for  a  2-D  function  in  the  APL 
environment  as  in  the  following  example: 
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2.  Restart  STATGRAF  and  select  the  section  “experimental  design.”  The  fol¬ 
lowing  will  be  displayed. 


X.  Full  and  Fractional  Factorials 

2.  Central  Compos its  Designs 

3.  Alias  Structure 

4.  Response  Surface  Plotting 


3.  Select  “Response  Surface  Plotting”  and  make  the  entry  in  the  menu  as  de¬ 
scribed  in  Figure  A.6. 


Response  Surface  Plotting 

X-axis  minimum:  0  *  Y-axis  minimum:  0 

X-axis  maximum:  IS  —*  Y-axis  maximum:  is 

— *  Humber  of  x-axis  intervals:  15  *  Number  of  y-axis  intervals:  15 

— •  Function  type:  Matrix  input 


Figure  A.6.  Surface  Plotting  Menu. 

4.  Press  the  F6  key  and  fill  in  the  entries  as  shown  in  the  example  of  Figure  A. 7. 


Response  Surface  Plotting 

X-axis  minimum:  0  Y-axis  minimum:  0 

X-axis  maximum:  16  Y-axis  maximum:  16 

Number  of  x-axis  intervals:  15  Number  of  y-axis  intervals:  15 

Function  type:  Matrix  input 

-*  Grid  matrix:  XX 

Figure  A. 7.  Selecting  the  Function  Definition. 
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Surface  Plotting  Options 


Top  title:  THE  MATRIX  XX  (EXAMPLE) 

C2  lines) 

X-Axis  Title:  X 
X-Axis  Title:  X 
2-Axis  Title:  2 

-*  Plot  type:  Surface 

For  surface  plots: 

X-axis  -esolution:  100 
X-axis  resolution:  100 
Lines  parallel  to:  Both  axes 

For  contour  plot: 

Type:  Lines 
Contour  levels: 
or 

Number  of  Divisions:  10  Minimum:  1 

Figure  A .8.  Surface  Plot  Entry. 


X-axis  skip  factor:  1 
X-axis  skip  factor:  l 
Hidden  line  removal:  None 


Maximum:  256 


TOE  MftTSlX  XX  (EXAMPLE 


300 
230 
200 
2  150 
100 
50 
0 


Figure  A.9.  Surface  2-D  Plot. 

5.  Press  the  F6  key  and  enter  the  values  for  the  surface  or  contour  plot  as  shown 
in  the  example  of  Figure  A. 8. 

6.  Press  the  F6  key  to  execute  the  plot,  the  F4  key  to  print  the  plot,  and  the  F5 
key  to  edit  the  plot  (See  Figue  A.9). 

7.  Press  the  F10  key  to  leave  this  screen  and  return  to  the  menu. 


aC9w3KSwS«!^ 


WARNING:  Since  the  number  of  points  in  3-D  plots  may  exceed  the  memory 
space  available,  do  not  select  hidden  or  partially  hidden  lines  in  complicated 
plots.  Select  the  entry  “points”  instead  of  entry  “lines”  when  drawing  compli¬ 
cated  contour  plots. 


APPENDIX  B 

DIGITAL  SIGNAL  PROCESSING  LIBRARY 


Listed  below  are  functions’  names  and  brief  discussions  of  the  various  computer 
algorithms  used  to  accomplish  this  thesis  research.  All  programs  were  written  in 
APL  by  the  author  unless  otherwise  noted.  Program  listings  and  samples  of  runs 
are  given  in  Appendix  C.  For  each  function  its  name  and  the  page  containing  its 
listing  in  Appendix  C  are  given.  Then  the  calling  sequence  is  given,  followed  by  a 
brief  description.  j 

A.  UTILITY  WORKSPACE 
1.  NORD  (Pg.  68^ 

Z«— N  NORD  S 

This  function  generates  a  normally  distributed  random  vector  using  the 
function  UNRD. 

The  function  inputs  are: 

a.  N:  Number  of  elements  in  the  vector 

b.  S:  Vector  with  two  arguments: 

*[S1]:  The  required  mean 
*[S2]:  The  required  standard  deviation 

The  function’s  source  is  the  software  package  GRAFSTAT. 


This  function  generates  a  uniformly  distributed  random  vector.  The  inputs 


are: 

a.  N:  Numbers  of  elements  in  the  vector 

b.  B:  Required  mean 


3.  PUTDATA  (Ps.  7(tt 

F  PUTDATA  ’A:\WORK\NAME. EXT > 

This  function  transfers  data  (vector,  array,  etc.)  in  free  format  from  an 
APL  workspace  into  an  IBM- PC  DOS  file.  The  fimction  uses  APL  NATIVE  files 
[7].  A  bell  sounds  when  the  data  transfer  has  been  completed.  The  function  inputs 


are: 

a.  F:  The  data  to  be  transferred 

b.  ’A:\WORK\NAME.EXT’:  The  designated  native  file  or  path:  “A:”  Stands  for 
the  drive  type,  “WORK”  stands  for  the  directory,  and  “EXT”  stands  for  extension 
of  the  DOS  file. 


WARNING:  The  user  is  not  permitted  to  designate  a  file  that  already  exists 
as  a  DOS  NATIVE  file. 


4.  GETDATA  (Pg.  72) 

F<— GETDATA  ’A:\WORK\NAME.EXT’ 

This  function  transfers  data  (vector,  array,  etc.)  from  an  IBM-PC  DOS 
file  to  an  APL  workspace.  The  function  uses  APL  NATIVE  files  [7].  A  bell  sounds 
when  the  data  transfer  has  been  completed. 

’  A :  \W0RK\NAME .  EXT  ’ :  The  native  file’s  name  or  path  has  to  be  transferred: 
“A :  ”  stands  for  the  drive  type,  “WORK”  stands  for  the  directory,  and  “EXT”  stands 
for  extension  of  the  DOS  file. 


B.  EC3400  WORKSPACE  [4] 


1.  DIGFRE 


TC«-FS  DIGFREQ  FC 

This  function  calculates  the  digital  frequency  using  the  equation: 


9  =  uT  =  2tt/  •  1  //, 


(B.  1] 


f  =  The  analog  frequency 
fa  =  The  sampling  frequency 
The  inputs  are: 

a.  FS:  The  sampling  frequency  [Hz] 

b.  FC:  The  analog  frequency  [Hz] 

The  output  is: 

TC:  The  digital  frequency  [rad.] 


2.  LPCOEFF  (Pg.  761 


h.LP't—n  LPCOEFF  TC 

This  function  calculates  the  21  +  1  lowpass  FIR  filter  coefficients  using  the 
Frequency  Transformation  and  the  equation: 


hLp(n)  =  V7r(n  ~  -0  sin([n  -  I}9C),  n  =  0, 1, 2, . . . ,  21 


(B.2) 


where  k  —  1  and  9C  =  The  digital  cutoff  frequency 

The  calculated  coefficients  axe  for  a  causal  nonrecursive  filter.  In  order  to 
generate  noncausal  coefficients  use  the  function  SHIFT.  The  inputs  are: 

a.  n:  Number  of  required  coefficients  (Use  an  odd  number.)  By  choosing  a  final 
number  of  coefficients  you  automatically  define  a  rectangular  window. 

b.  TC:  The  digital  cutoff  frequency  [rad.] 


3.  HPCOEFF  fPg.  77^ 
hHP<— n  HPCOEFF  TC 

This  function  calculates  the  21  +  1  highpass  FIR  filter  coefficients  using 
the  Frequency  Transformation  method  and  the  equation: 


hHP(n)  =  (-l)(n  7)hip(n),n  =  0,1,2,.. .  ,2/ 


(B.  3) 


The  calculated  coefficients  axe  for  a  causal  nonrecursive  (linear  phase) 
filter.  In  order  to  generate  noncausal  coefficients  (zero  phase  filter)  use  the  function 
SHIFT.  The  coefficients  are  normalized.  The  inputs  are: 

a.  n:  The  number  of  the  required  coefficients  (Use  an  odd  number.)  By  choosing 
a  final  number  of  coefficients  you  automatically  define  a  rectangular  window. 

b.  TC:  The  digital  cutoff  frequency  [rad.] 

4.  BPCOEFF  fPff.  78) 
hBP«—n  BPCOEFF  TUL 

This  function  calculates  the  21+1  bandpass  FIR  filter  coefficients  using 
the  Frequency  Transformation  method  and  the  equations: 

hBp(n )  =  [2  cos  [(n  -  I)90]\  hLP(n),  n  =  0, 1, 2, . . . ,  21  (BA) 

9C  =  (0U  ~  Bt)  /2 
9o  =  (9u  +  9e)/  2 

9U  —  The  upper  digital  frequency 
Of  =  The  lower  digital  frequency 
9C  =  The  digital  cutoff  frequency 

The  calculated  coefficients  are  for  a  causal  nonrecursive  filter.  In  order 
to  generate  noncausal  coefficients  use  the  function  SHIFT.  The  coefficients  are 
normalized.  The  inputs  to  the  function  are: 
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a.  n:  The  number  of  the  required  coefficients  (Use  an  odd  number.)  By  choosing 
a  final  number  of  coefficients  you  automatically  define  a  rectangular  window. 


b.  TUL:  Vector  of  two  arguments 

TUL  [1]:  The  upper  frequency  [rad.] 
TUL  [2]:  The  lower  frequency  [rad.] 


5.  BSCOEFF  fPg.  79) 
hBS<— n  BSCOEFF  TUL 

This  function  calculates  the  27  +  1  bandstop  FIR  filter  coefficients  using 
Frequency  Transformation  and  the  equation: 


hBs( 0)  =  1  -  Abp(0) 


hBs(n)  =  —hBp(n),  n  =  1, 2, . . . ,  27 


(B.  5) 


$c  =  (0u-0t)/2 

0o  =  (0u  +  9t)/ 2 

9U  =  The  upper  digital  frequency 
9t  —  The  lower  digital  frequency 
9C  =  The  digital  cutoff  frequency 

The  calculated  nonrecursive  coefficients  are  for  a  causal  nonrecursive  filter. 
In  order  to  generate  noncausal  coefficients  use  the  function  SHIFT.  The  coefficients 
are  normalized.  The  inputs  to  the  function  are: 

a.  n:  The  number  of  the  required  coefficients  (Use  an  odd  number.)  By  choosing 
a  final  number  of  coefficients  you  automatically  define  a  rectangular  window. 

b.  TUL:  Vector  of  two  arguments 


TUL  [1]:  The  upper  frequency  [rad.] 
TUL  [2]:  The  lower  frequency  [rad.] 


6.  Hamming  fPg.  80] 


W<— HAMMING  N 


This  function  generates  an  21  +  1  point  Hamming  window  using  the  equa- 


W(n)  =  0.54  +  0.46  cos  ^  ,  n  =  0, 1, 2, . . . ,  21 


(B.  6) 


The  window  starts  at  n  =  0  and  has  21+1  samples.  The  input  to  the  function  is 
N:  The  number  of  samples. 

7.  HANNING  (Ps.  81^ 

W HANNING  N 

This  function  generates  an  2 1  + 1  point  Hanning  (von  Hann)  window  using 


the  equation: 


W(n)  =  0.5  +  0.5  cos  ^  ,n  =  0, 1,2, . . . ,  21 


(5.7) 


The  window  starts  at  n  =  0  and  has  21+1  samples.  The  input  to  the 
function  is: 

N:  The  number  of  samples. 

8.  BARTLETT  (Pg.  82^ 

W+— BARTLETT  N 

This  function  generates  an  21+ 1  point  Bartlett  (triangular)  window  using 


the  equation: 


VT(n)  = 


1  -  }(n  -  I)  ,0  <  n  <  21 


(B.8) 


(  0  ,  otherwise 

The  window  starts  at  n  =  0  and  has  21+1  samples.  The  input  to  the 
function  is: 

N:  The  number  of  samples. 


9.  FREQRES  (Pg.  83^ 

HW<-Z  FREQRES  h 

This  function  generates  the  magnitude  and  phase  of  the  frequency  response 
from  the  unit  impulse  response.  The  functions:  FFT,  XMAGN,  XPHASE  are  used.  In 
order  to  achieve  good  resolution  the  coefficients  are  zero  padded  to  ’Z’  samples. 


WARNING:  The  ’Z’  needs  to  be  a  power  of  2. 


The  inputs  to  the  function  are: 

a.  Z:  The  toted  number  of  samples  including  the  zero  padding 

b.  h:  The  impulse  response  coefficients 
The  outputs  are: 

a.  HW:  The  frequency  response  magnitude  (Z/2  samples) 

b.  PW:  The  frequency  response  response  phase  (Z/2  samples)  [rad.] 


10.  IDEALF  (Pg.  841 
H+— TYPE  IDEALF  TC 

This  function  produces  the  frequency  response  of  various  ideal  filters.  The 
inputs  are: 

a.  ’TYPE’:  The  type  of  the  filter  to  be  generated: 

’LP’:  Lowpass 
’HP’:  Highpass 
’BP’:  Bandpass 
’BS’:  Bandstop 
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b.  TC:  A  vector  containing  the  following: 

1.  In  the  case  of  LP  or  HP: 

TC  [1]  =  The  number  of  samples 
TC  [2]  =  The  cutoff  digital  frequency 

2.  In  a  case  of  BP  or  BS: 

TC  [1]  =  The  number  of  samples 
TC  [2]  =  The  upper  digital  frequency 
TC  [3]  =  The  lower  digital  frequency 

The  frequencies  should  be  normalized  to  lie  in  the  interval  0  to  ir.  [rad.] 
11.  FCOEFF  (Pg.  861 
h<-n  FCOEFF  HW 

This  function  calculates  the  21+ 1  point  filter  coefficients  using  IFFT  (The 
sampling  method.) 

The  coefficients  axe  for  a  causal  nonrecursive  filter.  The  inputs  to  the 
function  are: 

a.  n:  The  number  of  the  required  coefficients  (Use  an  odd  number.)  By  choosing 
a  final  number  of  coefficients  you  automatically  define  a  rectangular  window. 

b.  HW:  The  ideal  frequency  response. 

C.  EC3410  WORKSPACE  [5,  8] 

1.  SACF  fPg.  881 
R<— LAG  SACF  DAT 

This  function  computes  the  sample  autocorrelation  function  (unbiased) 
using  the  following  equation: 

AT  — 1  — JJt| 

=  1/JV  £  r»+l*l  +y».  0<k<N-l  (B.  9) 

n=  0 
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The  inputs  are: 

a.  DAT:  The  data  block  (one-dimensional  vector) 

b.  LAG:  Number  of  lags  to  be  calculated  (the  limit  is  300  lags) 


2.  MEAN  (Pg.  90^ 
M<— MEAN  F 


This  function  calculates  the  sample  mean  of  a  data  set  using  the  following 


equation: 


N 


m  =  i/w£r<‘> 


(BIO) 


i=l 


The  input  is: 

F:  The  data  block 


3.  COY  (Pg.  91^ 
K<— COV  R 


This  function  generates  a  Toeplitz  covariance  matrix  using  covariance  func¬ 
tion.  The  input  is: 

R:  The  covariance  function  values  (vector) 

4.  LPFP  (Pg.  92) 
a«— P  LPFP  DATA 

This  function  generates  the  p th  order  Linear  Predictive  Filter  parameters 
and  calculates  the  Prediction  Error  Variance  by  solving  the  Normal  Equation: 


R(i)  R( 1)  ... 

T— 1 

1  ... 
a. 

'  1  ' 

a\ 

_ 

V o  •• 

1 _ 

i — 
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1  " 
»— * 

f?(0)  . 

a 

i  ‘ 
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The  inputs  are: 
a.  P:  The  coefficients  order 


b.  DATA:  The  data  block 


The  outputs  are: 

a.  a:  The  filter  coefficients 

b.  PRER:  The  estimate  prediction  error  variance 
The  following  functions  are  used: 

a.  SACF:  Calculates  the  sample  autocorrelation  of  the  data  block 

b.  COV:  Generates  the  Toeplitz  matrix  of  the  covariance  function 

5.  CC  (Pg.  931 
Y<— x  CC  h 

This  function  performs  one-dimensional  circular  convolution  between  two 
inputs  ( h(n )  @r(n)).  x  and  h  are  the  arguments  to  be  convolved. 


WARNING:  x  and  h  need  to  be  the  same  size  and  one-dimensional, 
otherwise  an  error  message  will  be  generated  and  the  convolution  will  not 
_ start. _ 

6.  FCC  (Pg.  941 
Y<— x  FCC  h 

This  function  does  one- dimensional  fast  circular  convolution  between  two 
arguments.  The  fast  operation  is  established  by  a  shifted  matrix  multiplication. 
x  and  h  are  the  arguments  to  be  convolved.  Due  to  memory  size  this  function 
is  limited  to  arguments  of  less  than  80  numbers.  If  the  argument  is  too  long  a 
printout  message  advises  the  user  to  use  regular  circular  convolution  function. 


WARNING:  x  and  h  need  to  be  the  same  size  and  one-dimensional, 
otherwise  an  error  message  will  be  generated  and  the  convolution  will  not 
start. 


w 

1 


-v r>'7> ' >’ 'j.'1^ T x ■  - 1  j.  •,■ v' v*  vr ■.■> v* 7T '.-» vV '.v '.m .'J '.m '.v.v 


7.  LCV  fPg.  95^ 

Y«— xl  LCV  hi 

This  function  does  one  dimensional  linear  convolution  between  two  argu¬ 
ments  by  using  the  following  equation: 

OO 

X\ (n)  *  hi(n)  =  ^2  xi(m)hi(n  —  m)  (J5.12) 

771=  —  OO 

xl  and  hi  are  the  arguments  to  be  convolved. 

The  function  uses  the  CC  function  (Circular  Convolution)  and  the  linear 
convolution  is  produced  by  zero  padding  the  arguments. 

8.  FLCV  fPg.  96) 

Y<— xl  FLCV  hi 

This  function  does  one-dimensional  linear  convolution  (same  as  LCV  func¬ 
tion)  by  using  the  function  FCC  in  order  to  achieve  a  faster  computation  time. 

9.  FFT  (Per.  971 

Z«-FFT  X 

This  function  calculates  the  one-dimensional  Discrete  Fourier  Transform 
according  to  the  definition 

N- 1 

X(k)=  J2  x(,n)e-^/N)nk,k  =  0,1,...,JV  — 1  (J3.13) 

n=0 

The  input  to  the  function: 

X:  Can  be  either  real  vector  of  length  N  or  a  two  by  N  matrix  (the  1st  row  is  the 
real  part  and  the  2nd  row  is  the  imaginary  part).  N  must  be  a  power  of  two. 

The  output  is: 


u 


[y: 


K: 


& 


«• 

ft 


I 


Z:  A  complex  two  by  N  matrix  (real  and  imaginary) 

The  function  was  written  by  Professor  Paul  Penfield,  Jr.  of  M.I.T. 


56 


10.  IFFT  (Pg.  981 

Z«— IFFT  X 

This  function  calculates  the  one- dimensional  Inverse  Discrete  Fourier  Trans¬ 
form  according  to  the  definition: 


x(n)  =  1/jV 


X(k)e-j(2*/N)nk 

.k=0 


,  n  =  0, 1, . . .  ,N  —  1 


(5.14) 


The  input  to  the  function  is: 

X:  Can  be  either  a  real  vector  of  length  N  or  a  two  by  N  matrix  (the  1st  row  is 
the  real  part  and  the  2nd  row  is  the  imaginary  part).  N  must  be  a  power  of 
two. 

The  output  is: 

Z:  A  complex  two  by  N  matrix  (real  and  imaginary) 


The  function  was  written  by  Professor  Paul  Penfield,  Jr.  of  M.I.T. 

11.  PSE  (Pg.  100) 

INW<— ZP  PSE  DATA 


This  function  computes  the  power  spectrum  estimate  (Periodogram)  using 
the  equation: 


Sz  ( eJ 6)  =  IN(6)  =  l/N 


X  (z]6) 


(5.15) 


The  function  uses  the  FFT  APL  function  and  the  inputs  are: 

a.  DATA:  The  data  block  (one  dimensional  vector) 

b.  ZP:  The  length  of  the  zero  padding  needed  for  a  smooth  plot 


WARNING:  The  DATA  plus  ZP  need  to  be  in  length  which  is  a  power  of  two 
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12.  PSEB  (Pg.  100) 

INBW«— SEG  PSEB  DATA 
This  function  computes  the  power  spectrum  estimate  with  rectangular 
nonoverlapping  window  using  the  equation: 


IN'(i){$)  =  1  /N  X{i)  ( eje ) 


M 

B(0)  =  1  /M^2lN'(i)(0) 

i=i 


(B.  16) 


N:  The  length  of  the  data  set 
M:  The  length  of  a  segment 

Each  segment  has  N'  Points:  N  =  N'  x  M. 

This  function  uses  the  FFT  and  the  PSE  APL  functions  and  the  inputs  are: 

a.  SEG:  The  length  of  the  segments 

b.  DATA:  The  data  block  (one-dimensional  vector) 

The  function  automatically  zero  pads  the  segment  up  to  256  points. 


WARNING:  In  case  of  error  due  to  memory  size,  use  fewer  segments 


14.  XTIME  (Pg.  1021 


XZ«— XL  XTIME  XR 


This  function  multiplies  two  complex  numbers 


XZ  =  XL  x  XR 


(B.  18) 


15.  XMAGN  (Pg.  1031 


Z*— XMAGN  XR 


This  function  computes  the  complex  magnitude  of  a  number,  the  result  is 


a  real  number 


Z  =  |XR| 


(£.19) 


16.  XPHAS  (Pg.  1031 


Z<— XPHAS  XR 


This  function  computes  the  phase  of  the  complex  number  XR,  and  returns 


a  real  number. 


Z  XR 


(B.  20) 


17.  XCONJ  (Pg.  1041 


XZ< — XCON J  XR 


This  function  returns  the  complex  congugate  of  the  input. 

18.  XEPO  (Pg.  104) 


XZ<~XEP0  XR 


This  function  computes  the  complex  exponential 


XZ  =  e" 


(£-21) 


19.  PI  (Pg.  1051 


Y<— PI 


k  *M»>  M  *  *  *  •  — f  *  *  *  •  J J*"*  ***  »•  *»  *  ■  “  »  ’  •*  "  •»  j*»  ■  «  ^ ^  M*_  MV  A.  JV  «/■  « v  mV  ^ 


This  function  returns  the  7r  value. 


20.  COS  (Pet.  105) 

Y*— COS  X 

This  function  returns  the  cosine  of  the  input. 

21.  SIN  (Pa.  1051 

Y«— SIN  X 

This  function  returns  the  sine  of  the  input. 


D.  EC4440  WORKSPACE  [9] 
1.  CC2D  fPg.  1071 
Y<— X  CC2D  H 


This  function  calculates  two  dimensional  circular  convolution  using  the 


equation: 


Y  =  h  *  *  x 

Nx~  1  AT2-1 

Y{nun2)~Y^  ^2  h(nun2)  x[((m  -m1))Nl((n2 -m2))N2]  {B.  22) 

nl  =  0  n2=0 

Where  subscript  Ni  and  N 2  denotes  arithmetic  modulo  N\  and  N 2.  x  and  h  are 
the  arguments  to  be  convolved. 


WARNING:  X  and  H  need  to  be  two-dimensional  with  the  same  size,  other¬ 
wise  an  error  message  will  be  generated  and  the  convolution  will  not  start. 


2.  FFT2D  fPg.  1091 

Y<— FFT2D  X 

This  function  calculates  the  two-dimensional  DFT  using  the  one-  dimen¬ 
sional  FFT  function.  If  the  input  is  real,  the  function  converts  it  to  a  complex  array. 
If  the  output  array  is  real  valued  the  second  plane  is  eliminated  automatically. 


.A- i 


,  V.V.  / 


,  %  n 


60 


WARNING:  The  input  needs  to  be  a  length  which  is  a  power  of  two. 


3.  IFFT2D  (Pg.  llCh 

Y<— IFFT2D  X 

This  function  calculates  the  two-dimensional  Inverse  DFT  using  one-dimensional 
IFFT  function  (the  same  explanation  as  FFT2D). 

4.  SHIFT  (Pg.  1141 

Y<— SHIFT  X 

This  function  shifts  periodic  data  to  an  interval  starting  at  the  origin. 

5.  UNSHIFT  (Pg.  1141 

Y<— UNSHIFT  X 

This  function  shifts  periodic  data  to  an  interval  centered  at  the  origin. 

6.  SHIFT2D  (Pg.  115) 

Y<— SHIFT2D  X 

This  function  shifts  two-dimensional  periodic  data  to  an  interval  starting 
at  the  origin 

7.  UNSHIFT2D  (Pg.  1151 

Y<— UNSHIFT2D  X 

This  function  shifts  two-dimensional  periodic  data  to  an  interval  centered 

at  the  origin.  Each  of  the  following  functions  generates  a  two-dimensional  window 

by  using  a  one-dimensional  window.  Two  types  of  supports  are  used: 

a.  A  rectangular  region  of  support  which  is  achieved  by  computing  the  outer 
product  of  two  1-D  windows: 

WR(ni,n2)  =  W,(n,)W2(n2)  (B.  23) 
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b. 


A  circular  region  of  support  which  is  formed  by  sampling  a  circularly  rotated, 
1-D,  continuous  function: 


Wc{nun2)  =  w(^^+n\Sj 


(B.  24) 


8.  ILPFILT  fPg.  117} 

W<-R  ILPFILT  P 

This  function  generates  a  circular  based  rectangular  window  using  the  1-D 
equation: 

Vr(n)  =  It  ’  “  (B.  25) 


i\  <  R 
otherwise 


-(.«)= {1  ,|n! 

;  \  0  ,  otl 

This  can  also  be  used  for  defining  the  frequency  response  of  an  ideal  circular  sym¬ 
metrical  lowpass  filter.  The  inputs  axe: 

a.  R:  The  window  base  radius 

b.  P:  The  dimension  of  the  square  matrix  in  which  this  window  is  defined. 


9-  HAMMINGR  fPq.  1191 

WR<— II  HAMMINGR  N 

This  function  generates  a  2-D  Hamming  rectangular  based  window  using 
the  1-D  equation: 


W(n)  =  {  0,54  ~  0,46  cos  (  W-0  ,  M  <  I 

(  0  ,  otherwise 


(5.26) 


The  inputs  are: 

a.  II:  The  window  base  dimension  (one-side) 

b.  N:  The  dimension  of  the  square  matrix  in  which  this  window  is  defined. 


10.  HAMMINGC  (Per.  121) 

WC«-II  HAMMINGC  N 

This  function  generates  a  2-D  Hamming  window  with  circular  base  using 
equation  (B.26).  The  inputs  are: 

a.  II:  The  window  base  radius 

b.  N:  The  dimension  of  the  square  matrix  in  which  the  window  is  defined. 


11.  HANNINGR  fPg.  123) 

WR+— II  HANNINGR  N 

This  function  generates  a  2-D  Hanning  window  with  rectangular  base  using 
the  1-D  equation: 


0.5  [1  +  cos  (7 rn/I)\ 
0 


,  |n|  <  / 

,  otherwise 


(B.  27) 


The  inputs  are: 

a.  II:  The  window  base  dimension  (one  side) 

b.  N:  The  dimension  of  the  square  matrix  in  which  the  window  is  defined. 

12.  HANNINGC  (Pg.  125^ 

WC<-II  HANNINGC  N 

This  function  generates  a  2-D  Hanning  window  with  circular  base  using 
equation  (B.27).  The  inputs  are: 

a.  II:  The  window  base  radius 

b.  N:  The  dimension  of  the  square  matrix  which  this  window  is  defined. 


13.  BARTLETTR  (Ps.  1271 
WR<— II  BARTLETTR  N 
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This  function  generates  2-D  Bartlett  (triangular)  window  with  rectangular 


base  using  the  1-D  equation: 


1  —  \n  ,  0  <  n  <  I 
W(n )  =  <  1  +  jn  ,  — 1  <  n  <  0 
k  0  ,  otherwise 


(5.28) 


The  inputs  are: 


a.  II:  The  window  base  radius 

b.  N:  The  dimension  of  the  square  matrix  in  which  the  window  is  defined. 


14.  BARTLETTC  fPg.  129) 


WC«-II  BARTLETTC  N 

This  function  generates  a  2-D  Bartlett  window  with  rectangular  base  using 
equation  (B.28).  The  inputs  axe: 

a.  II:  The  window  base  radius 

b.  N:  The  dimension  of  the  square  matrix  in  which  the  window  is  defined. 


15.  RECTANR  (Pg.  131^ 


WR«— II  RECTANR  N 

This  function  generates  2-D  rectangular  window  with  rectangular  base 
using  equation  (B.25).  The  inputs  are: 

a.  II:  The  window  base  dimension  (one  side) 

b.  N:  The  dimension  of  the  square  matrix  in  which  the  window  is  defined. 


16.  PROTFILT  fPg.  133) 


HPROT*— PROTFILT  SIZE 

This  function  produces  a  frequency  response  1-D  lowpass  prototype  filter 
in  order  to  design  a  2-D  FIR  filter  using  a  transformation.  The  function  SHIFT 
has  been  used.  The  input  is: 


v.v.v 


W5S 


SIZE:  The  size  of  the  prototype  filter  (number  of  samples). 


17.  COEFF  (Per.  1351 

A<— COEFF  HP 

This  function  finds  the  coefficients  a(n)  in  the  representation: 


N 


h(u)  =  y,  a(n)  cos(am) 


n=0 


where: 


,  ,  A  f  h(o)  ,  n  -  0 
a(n)  —  (  2h(n)  ,n  >  0 


(B.  29) 


The  function  IFFT  has  been  used.  The  input  is: 

HP:  The  filter’s  1-D  prototype. 

18.  TRANSFNC  (Pg.  1371 

Fww*— TRANSFNC  SIZE 

This  function  produces  a  transformation  matrix  to  design  a  FIR  filter 
using  McCLELLAN  transformation.  In  this  example  the  following  function  has 
been  used: 


F(uji  ,u>2)  —  1/2  (  —  1  +  coswi  +  cosu;2  +  coslji  cosu^)  (R.30) 

In  order  to  change  F(ui j ,  u>2 )  to  a  different  transformation  function,  replace 
lines  10,  13,  and  14  in  the  program  with  the  new  choices.  (See  example  in  Appendix 
C  (Pg.  137).  The  input  is: 

SIZE:  The  number  of  samples  of  u)\  or  u>2. 

19.  CHEB  fPg.  1401 

Y<-N  CHEB  X 
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This  function  calculates  the  Chebyshev  Polynomials  using  the  following 


recursion: 

T„(x)  =  2xTn_i(x)  -  Tn_2(x)  (B. 31) 

The  inputs  are: 

a.  N:  The  order  of  the  desired  polynomial 

b.  X:  The  argument. 


WARNING:  The  combination  of  high  order  of  Chebyshev  polynomial  and  a 
large  array  to  be  evaluated  will  result  in  a  long  period  of  calculating  time. 


20.  MCCLEL  (Pg.  141>) 

Hww*— a  MCCLEL  F 

This  function  computes  the  2-D  frequency  response  using  McCLELLAN 
transformation  according  to  this  equation: 

N 

H(u i,w2)=  Y  a(n)Tn  [F(u;i,u;2)]  {B. 32) 

n=0 

In  order  to  indicate  the  calculating  status,  the  function  prompts  the  present  cal¬ 
culating  order  of  the  polynomial.  The  inputs  are: 

a.  a:  The  impulse  response  coefficients  of  the  prototype  filter 

b.  F:  The  transformation  matrix  (previously  calculated  using  the  function  TRANSFNC) 
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APPENDIX  C 

FUNCTION  LISTINGS  AND  EXAMPLES  OF  USE 


3 LOAD  0  UTILITY 

0  UTILITY  SAVED  12/10/1987  16:37:45 

I  NAVAL  POSTGRADUATE  SCHOOL 


MONTEREY. CALIFORNIA  I 


VERSION  1.0 


UTILITIES  WORKSPACE 

THESIS  BY 
Y.KATZIR.I.A.F. 

ADVISORsPROF.  C.W.  THERRIEN. 

SEPTEMBER  1987 


GETDATA  INFORMATION  NORD  PUTDATA  UNRD 


W.V.V.V.N ' 


,  </v' 


VNORDCtnv 

to:  z«-n  nord  P:s:i:T 
Cl]  ft 

C2 J  ft  A  NORMAL  DISTRIBUTION  RANDOM  VECTOR  GENERATOR 
C  3  ]  ft 

C4  ]  P«*2TP,1 

C5]  Z«-CN)pO 
C6]  I<-1 

C7]  F10 : T«-2  ONRD  0.5  0.5 

C8]  T«-C2xT)-l 
C9]  S«-(TC1]*2)+TC2]*2 

CIO]  -»F10xT.sai 

Cll]  ZCI]*-PCl]+PC2]xCTCl]xCC_2xftS]^S]*0.5] 

C 12  ]  -*F10xi.  CI*-I*1)  aN 

C 13 ]  *0 

C 14  ]  ft 

C15]  ft  SOURCE:  GRAFSTAT 

CIS]  ft  N«- NUMBER  OF  ELEMENTS  IN  THE  VECTOR 

C 17  ]  ft  P*-MEAN 

C 18  ]  ft  S«-SIGMA  CNOT  SIGMA*2] 

C 19]  ft 


VUNRDC0] V 
CO]  R«-N  UNRD  B 
Cl]  ft 

C2]  ft  A  UNIFORM  DISTRIBUTION  RANDOM  VECTOR  GENERATOR 
C3]  ft 

C  4 ]  R«-CBC1]-BC2])  +  CCN?10000]x2xBC2])+10000 

C5]  -*0 

C  6]  ft 

C7]  ft  SOURCE:  GRAFSTAT 

C  8  ]  ft  N«- NUMBER  OF  ELEMENTS  IN  THE  VECTOR 
C9]  ft  B«-THE  LIMITS  C2  NUMBERS) 

CIO]  A 
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GAUSEw-100  NORD  0  1 


GAUSD 

1.558274527  "2.444846636  "1.09846582  1.122362631  0.580943715 
"0.2714227577  0.4143999987  '0.9733729641  "1.021759908 
0.3179468935  1.515647616  0.7495965254  "0.5068186108  0.8851652886 
"0.2479112202  "0.7262416263  "0.4452313973  "0.6121927891 
"0.2086518602  0.5618385258  "1.064013245  0.3513946699  1.132830543 
0.1510880125  0.7028398161  "0.05231570544  2.018906718 
0.9240866379  "1.814329439  0.0351494428  "1.806312062  1.028100377 
0.3944470661  0.6391243468  0.3738735133  1.75599265  “0.3197284973 
"0.1366870369  0.6159939661  0.9774365509  "1.115427909 
"0.5496175242  0.03996605725  "2.484219992  1.158618634 
"1.026625814  1.153535117  "0.7858857486  0.6349969194  0.8198888762 
"0.1760331562  0.5627148978  "0.127309154  0.5538766995  "1.09726933 
"0.7311467753  1.404690913  "0.6203249336  0.2370824053 
"1.585608223  "0.4013458241  "0.7707101307  "0.2637646461 
•  0.9762420248  0.9775308888  1.169573915  0.1594895548  0.5001385714 

"1.05500756  "0.4509335123  1.271139261  0.8989835792  0.438888127 
"1.247429466  0.3241026886  0.3904619568  "0.4052863418 
0.2929672888  2.565225342  "0.4576884119  "1.611439324  "2.669649901 
"0.7594268721  "0.6750622611  "1.171739437  2.032501263 
0.9683033719  0.6699827805  0.4199990777  "2.871065022  1.686296301 
0.02722627758  "0.9020480504  "2.053207747  0.08940644014  2.0866291 
0.3651320115  0.8457720835  "0.1842246182  1.03036774 


UNIFIX-100  UNRD  0  100 


UNIFD 

3.22  "5.78  18.74  "10.04  43.86  91.44  "63.84  "19.84  "18.58  94.66  "51.14 
"46.78  "80.44  "23.98  "67.32  69.64  71.06  "82.22  "14.14  "37.6 
"66.84  70.92  "93.12  "30.38  18.1  92.06  50.94  "90.96  82.12  24.6 

29.6  83.32  "1.1  34.96  "92.3  61.72  "81.84  91.02  "77.52  3.34  "18.2 
26.18  "16.3  "93.02  "75.96  "14.68  66.88  44.58  "75.88  "2.82  "10.68 
"87.08  "34.64  81.86  "61.82  "68.26  48.36  35.22  6.46  "76.02  "84.16 

26.7  74.7  "62.86  "31.94  0.42  "42.96  96.16  33.74  "52.68  39.3 
"51.16  "8.2  "77.86  "60.26  "3.18  "13.76  16.38  64.04  "5.46  "24.86 
"31.52  "58.7  "35.88  "16.94  "47.26  73.98  "87.34  "94.42  '62.04 
"80.84  61.58  "89.98  57.78  "99.9  43.1  60.76  "32.04  "1.38  "89.56 
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YPUTDATACD]? 

CO]  7  PUTDATA  NAME;OIO;TIE;FA:NREC;C:FB;FC;B£LL;I 
[II  a 

C2 1  A  FUNCTION  TO  WRITE  DATA  CF3  C VECTOR. AREA*. ETC. JIN  FREE  FORMAT 
C 3  3  A  INTO  AN  IBM  P.C.  D.O.S.  FILE  t ’NAME’ ] 

1 4  3  A 

CS  3  aiO«-I«-l 
C  6  3  QNUNTIE  "1 

C7]  TI£«-"1 

C8]  A  CREATING  A  TIED  FILE 
C91  NAME  ONCREATE  TIE 
CIO 3  A  WRITING  OUT  THE  SHAPE  OF  F 
till  FA*-".  10  0  «pF 

t 123  A  APPENDING  THE  SHAPE  OF  F  TO  THE  FILE 

C13]  < OTCLF . OTCNL . FA ) ONAPPEND  TIE 

C14]  A  RESHAPING  F 

C1S3  NREC*-rcx/pF)+4 

1 163  C«-1+4I  Cx/CpFJ)-l 

C17]  F*-(NREC,4]pF 

d8]  A  WRITING  OUT  MOST  OF  THE  ARRAY  F 
C19]  FB«-»C*1  0  IF) 

1201  LOOP: FBCI: CFBCI; ]* ' “' ]  /ipFBCI: 11+’-' 

C21]  -*CNREC>I*-I*l]/LOOP 

C22  3  C OTCLF . OTCNL . FB ) UNAPPEND  TIE 

C23 ]  A  WRITING  OUT  THE  REMAINDER  OF  F  AS  THE  LAST  RECORD 
C24  ]  FC«-'  *  ,  4CCTFCNREC;  ]  ] 

C25]  FCt CFC3 ]/\.pFC]«*'-' 

[26]  C OTCLF. OTCNL. FC1 ONAPPEND  TIE 

C27]  A  UNTYING  THE  FILE 

C28]  DNUNTIE  ~1 

C29  3  BELL«-QTCBEL 

C30  3  BELL 

C  3 1  ]  -*0 

[32  3  A 

C33I  A  THE  FUNCTION  OSES  NATIVES  APL  COMMANDS. 

134  3  A 

C3S3  A  Y.KATZIR  .I.A.F. .  6.1.87. 

C36J  A  ADVISOR:  C.W.  THEHRIEN. 

C 37 ]  A 
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/.V.V.'A'  '^'-k  v^^^->r^>,1'-x",-»  '  >■  ^  -,  w^i 
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SAOSD  PtJTDATA  'oAas.nAT' 

)ory 

C:\>«PB  SACS.DAT 


too 

1.558274327 
0.380943713 
-1.021739908 
-0.3068138108 
-0.4432313973 
-1.064013245 
0.7028398161 
-1.314329439  ‘ 
0.3944470661 
-0.3197284973 
-1.113427909 
1.158613634 
0.6349969194 
-0.127309134 
1.404690913 
-0.4013438241 
0.9773308888 
-1.05300736 
0.438888127 
-0.4052863418 
-1.611439324 
-1.171739437 
0.4199990777 
-0.9020480504 
0.3651320115  0. 


-2.444846636 

-0.2714227377 

0.3179468933 

0.3851652886 

-0.6121927391 

0.3313946699 

-0.05231570344 

0.0351494428 

0.6391243468 

-0.1366870369 

-0.5496173242 

-1.026623814 

0.3198888762 

0.5338766993 

-0.6203249336 

-0.7707101307 

1.169573915 

-0.4509333123 

-1.247429466 

0.2929672388 

-2.669649901 

2.032501263 

-2.871065022 

-2.053207747 


-1.09846582 
0.4143999987 
1.515647616 
-0.2479112202 
-0.2086518602 
1.132330343 
2.018906713 
-1.806312062 
0.8738733133 
0.6139939661 
0.03996603725 
1.133535117 
-0.1760331562 
-1.09726933 
0.2370824053 
-0.2637646461 
0.1394893548 
1.271139261 
0. 3241026886 
2.565225342 
-0.7594268721 
0.9683033719 
1.686296301 
0.08940644014 


1.122362631 
-0.9733729641 
0.7495965254 
-0.7262416263 
0.5613385258 
0.1510880125 
0.9240866379 
1.028100377 
1.75599263 
0.9774365509 
-2.484219992 
-0.7858857486 
0.5627148978 
-0.7311467733 
-1.385608223 
0.9762420248 
0.5001385714 
0.8989835792 
0.3904619568 
-0.4376884119 
-0. 6730622611 
0.6699827305 
0. 02722627758 
2.0866291 


8457720835  -0.1842246182  1.03036774 
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vn’av 


vgetdataco]? 

CO]  F«-GETDATA  NAME : OIO; DATA ;M;SHAP; I : BELL 

ci:  a 

C  2  ]  A  FUNCTION  TO  READ  FREE  FORMATED  DATA  CF]  CVECTOR.ARRAY.ECT. ) 

C3]  A  FROM  IBM  P.C.  D.O.S.  FILE  C ’ NAME ' 3 

u:  a 

C5]  0IO«-I«-0 

C 6 ]  ONUNTIE  "1 

C7]  A  TYING  THE  FILE 
C8]  NAME  ONTIE  ~1 

C9]  A  FINDING  THE  SI2E  OF  THE  FILE. 

CIO]  M«-0NSIZE  "1 

CllJ  A  READING  THE  DATA  INTO  THE  WORK  SPACE 

C12]  FH3NREAD  "1.82.M 

C 13  ]  FC  CF=0AVC45 ]  )  /l.pF]<-QAVC  253 1 

C 14  3  A  EXECUTING  THE  DATA  AS  AN  APL  EXPRESSION 

CIS]  F«-*('*CF=QAVC  10]  0  v (F^OAVC 13  ] ) WCF=QAVC26 ) ) ) /  F 

CIS]  ONUNTIE  “1 

C 17  ]  LOOP:  I+-I+1 

C 18 ]  -K 1210) /ERROR 

C19  ]  DATA«-I1F 

C  20  ]  SHAP*-ITF 

C21]  •♦(  Cx/SHAP)*CpDATA)  )/LOOP 

C  22  ]  F«-SHAPp DATA 

C23  ]  BELL«-0TC8EL 

C24 ]  BELL 

C25]  -*0 

C26]  ERROR: 'FILE  INCORECTLY  FORMATED' 

C27  ]  -*0 

'28]  A 

129]  A  DATA  HAVE  BEEN  CREATED  BY  APL, PL/ I. FORTRAN. EDITOR  ETC. 

C  30  ]  A 

C  31]  A  NUMBERS  NEED  NOT  HAVE  DECIMAL  POINTS  AND  ARE  DELIMITED  BY  SPACE. 
C 32 ]  A  FIRST  RECORD  IS  SHAPE  OF  DATA. 

C 33 ]  a 

C 34 ]  A  Y.KATZIR.I.A.F. ,  6.1.87. 

C 35 ]  A  ADVISOR:  C.W.  THERRIEN. 

C  36  ]  A 
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v  Irrrrrrrf  ■.*  <'  swjvivw.yv*  v*tw7  '^v.vw.  t  r.  a  k  w  wwv.-  ^ *.-  *.-  v  «-.  *.-.*.  ~ . 


GAUSCM-GETDATA 


’ GAOS . DAT ' 


f 

i 

t 

GAUSD 

1.558274527  "2.444846636  “1.09846582  1.122362631  0.580943715  "0.27142275 
0.4143999987  "0.9733729641  "1.021759908  0.3179468935  1.515647616 
0.7495965254  "0.5068186108  0.8851652886  "0.2479112202  "0.726241626 
[  "0.4452313973  "0.6121927891  *0.2086518602  0.5618385258  "1.06401324 

,  0.3513946699  1.132830543  0.1510880125  0.7028398161  "0.05231570544 

|  2.018906718  0.9240866379  “1.814329439  0.Q351494428  "1.806312062 

i  1.028100377  0.3944470661  0.6391243468  0.8738735133  1.75599265 

«  "0.3197284973  “0.1366870369  0.6159939661  0.9774365509  "1.115427909 

I  “0.5496175242  0.03996605725  “2.484219992  1.158618634  “1.026625814 

1.153535117  "0.7858857486  0.6349969194  0.8198888762  "0.1760331562 
i  0.5627148978  "0.127309154  0.5538766995  "1.09726933  "0.7311467753 

.  1.404690913  "0.6203249336  0.2370824053  "1.585608223  "0.4013458241 

>  "0.7707101307  “0.2637646461  0.9762420248  0.9775308888  1.169573915 

■  0.1594895548  0.5001385714  "1.05500756  "0.4509335123  1.271139261 

I'  0.8989835792  0.438888127  “1.247429466  0.3241026886  0.3904619568 

1*0.4052863418  0.2929672888  2.565225342  "0.4576884119  “1.611439324 

“2-669649901  *0.7594268721  “0.6750622611  "1.171739437  2.032501263 
0.9683033719  0.6699827805  0.4199990777  "2.871065022  1.686296301 
0.02722627758  “0.9020480504  “2.053207747  0.08940644014  2-0866291 
0.3651320115  0.8457720835  "0.1842246182  1.03036774 
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) LOAD  0  EC3400 

0  EC3400  SAVED  12/12/1987  21:23:23 


NAVAL  POSTGRADUATE  SCHOOL  MONTEREY .  CALIFORNIA 


ECE  -3400  WORKSPACE 

THESIS  BY 
Y.KATZIR.I.A.F. 

ADVISOR:  PROF.  C.W.  THERRIEN. 

VERSION  1.0  NOVEMBER  1987 


)FNS 

BARTLETT  BPCOEFF  BSCOEFF  CC  COS  DIGFREQ  FCC  FCOEFF  FFT 

FLCV  FREQRES  HAMMING  HANNING  HPCOEFF  IDEALF  I  FFT  INFORMATION  LCV 

LPCOEFF  MEAN  PI  SHIFT  SIN  UNSHIFT  XCONJ  XEXPO  XMAGN  XPHAS 

XSQRT  XTIME 


7DIGFREQEC33Y 
£03  TC«-FS  DIGFREQ  FC 
£13  A 

£23  A  CALCULATING  THE  DIGITAL  FREQUENCY 
£33  A 

£43  TC«-(2*PI*FC)+FS 
£53  -*0 

£63  A 

£7  3  A  FS ♦■SAMPLING  FREQUENCY  [Hz 3 

C8  3  A  FCt-THE  CORRESPONDING  ANALOG  FREQUNCY  £Hz3 

£93  A  TC:T!£E  DIGITAL  FREQUENCY  ERAD.  3 

£103  A 

£113  A  Y.XATZIR.I.A.F.  .NOVEMBER  1987 
£121  A 


FC*10000 

7S-50000 

TC«-FS  DIGFREQ  FC 


VLPC0EFFCD3V 

CO]  hLP*n  LPCOEFF  TC:I;QIO 
Cl]  A 

C2]  A  CALCULATING  THE  '!»’  LGWPASS  FUl  FILTERS  COEFFICIENTS  USING  FOURIER 

C31  A  METHOD 

C4]  A 

C5]  010*1 

C6]  n*n-l 

C7]  n*lCn+2) 

C  3  ]  hLP*npO 

C9]  hLP*(l+(PIi«(  in)  J  )xSINC  C  in)*TC) 

CIO]  hLP*(TC+PI) ,hLP 
cil]  hLP*(ntC«liLP) )  ,hLP 
C12]  010*0 

C13)  -*0 

C  14  ]  A 

CIS]  A  THE  COEFFICIENTS  ARE  FOR  CAUSAL  NONRECURSIVE  FILTER. 

CIS]  A  THE  COEFFICIENTS  ARE  NORMALIZED. 

C 17 ]  A  rn-ODD  NUMBER  OF  REQUIRED  COEFFICIENTS  C RECTANGULAR  WINDOW) 

CIS]  A  TC*OIGITAL  CUTOFF  FREQUENCE  [RAD.  ] 

C191  A 

C201  A  E.XATZIR. I. A.F.  .NOVEMBER  1987 
C21]  A 


TC*0.4*PI 

n*21 


hLP*n  LPCOEFF  TC 
tlLP 

”1.559907461E*17  "0. 03363S7435  ”0.02338723209  0.02672826525  0.05045511524 
”1.3599 07 4 6 IE *17  ”0.07568267286  *0.06236595225  0.09354892838 
0.3027306915  0.4  0.3027306915  0.09354892838  ”0.06236595225 
”0.07368267286  ”1. 559907461£”17  0.05045511524  0.02672826525 
”0.02338723209  ”0.0336367435  ”1. 559907461E”17 
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v  V  V  V  V.  <■„  /.  ■''.I-’* A  .  A/,1'' 


w 


C03  h 
Cl]  a 
C2]  a 
C33  A 
C4]  A 
C5I 
C6] 

C7: 

C  3  ] 

C  9 1 
CIO] 
cii:  a 
C121  A 
C 13  ]  A 
C14]  A 
CIS]  A 


THPCOEFFCD]? 

hHP<-n  HPCOEFF  TC:aiO:N:hLP 
A 

A  CALCULATING  THE  'n'  HIGHPASS  FIR  FILTER  COEFFICIENTS  USING  FOURIER 
A  METHOD 
A 

aioo 

N«-n 

TC«-PI-CTC) 
hLP«-n  LFCOEFF  TC 
hHP«-CC-l]*lN)xilLP 
•*0 


C12]  A  THE  COEFFICIENTS  ARE  FOR  CAUSAL  NONRECURSIVE  FILTER. 

C131  A  THE  MAGNITUDE  IS  NORMALIZED. 

C14]  A  n«-ODD  NUMBER  OF  REQUIRED  COEFFICIENTS  (RECTANGULAR  WINDOW) 
CIS]  A  TO-DIGITAL  CUTOFF  FREQUENCY  CRAD.  ] 

CIS]  A 

C17]  A  Y. XATZIR. I. A. F.  .NOVEMBER  1987 


TO-0.4KPI 

n«-2l 


hHP«-n  HPCOEFF  TC 


"2. 339861191E"17  0.0336367435  0.02338723209  "0.02672826525  "0.05045511524 
*2. 33986119 IE "17  0.07568267286  0.06236595225  "0.09354892838 
"0.3027306915  0.6  "0.3027306915  "0.09354892838  0.06236595225 
0.07568267286  "2.339S61191E"17  "0.05045511524  "0.02672826525 
0.02338723209  0.0336367435  "2. 339861191E"17 


. 


fbpcoeffcd]* 

C01  hBP*-n  BPCOEFF  TUL:OIO:TC;TO;hLP 
Cl]  A 

C2]  A  CALCULATING  THE  'n'  BANDPASS  PIR  FILTER  COEFFICIENTS  USING  FOURIER 

C3  3  A  METHOD 

C4]  A 

CS]  QI0«-0 

C6J  TC4-<TULC03-TULCl])+2 
C7]  T04-(TULC01+TULCl])+2 

CS]  HLP4-n  LPCOEFF  TC 
C9]  n«-n-l 
CIO]  n*lCn+2) 

Cll]  hLP4-nihLP 

C12]  hBP4-C2xC0SCCvn-H)xT0))xhLP 
C13  ]  hBP4-CnTC<»hBP]  ]  ,hBP 
Cl*]  -*0 
CIS]  A 

C16]  A  THE  COEFFICIENTS  ARE  FOR  CAUSAL  NONRECURSIVE  FILTER. 

C17]  A  THE  COEFFICIENTS  ARE  NORMALIZED. 

CIS]  A  n«-ODD  NUMBER  OF  REQUIRED  COEFFICIENTS  C RECTANGULAR  WINDOW] 

C19 ]  A  TUL4-0IGITAL  UPPER  AND  LOWER  CUTOFF  FREQUENCIES  (TWO  ELEMENTS)  CRADJ 
C20]  A 

C21]  A  1. XATZIR, I. A.F. .NOVEMBER  1987 
C22]  A 


TUL4-2pO 

TULCl]4-0.6xpi 

TULC2]4-0.*xPI 

n*-2l 


hBP*n  8PC0EFF  TUL 


hap 

"7. 799 5373058*1 8  1. 20*3S3817E“17  0.0*677446*19  *3 . 152585957E*17  *0.1009102 
3. 89 56273 05E *17  0.1513653457  *3 . 1568404598*17  *0.1870978568 
1. 207227 699E*17  0.2  1. 207227699E*17  *0.1370978568  '3.1568404598*17 
0.1513653457  9 . 895627305E*17  *0.1009102305  *3 . 152585957E"17 
0.04677446419  1. 204383317E*17  *7.799537305E*18 


•  wsrs 


fbscoeffco:* 

hBS«-n  8SCOEFF  TUL:QIO;hflP 

A 

A  CALCULATING  THE  *n'  BANDSTOP  FIB  KILTER  COEFFICIENTS  USING  FOURIER 


csj  aio«-o 

C61  hflP<-n  BPCOEFF  TUL 
C7  ]  hBS*~ilBP 

C8]  n«-n-l 

C9]  hBSCn+23«-l-hBPCn+2] 

CIO!  -*0 

C113  A  THE  COEFFICIENTS  ABE  FOR  CAUSAL  NONRECURSIVE  FILTER. 

C12]  A  THE  COEFFICIENTS  ARE  NORMALIZED. 

C131  A  n«-ODO  NUMBER  OF  REQUIRED  COEFFICIENTS  (RECTANGULAR  WINDOW) 

C14  ]  A  TUL«-OIGITAL  UPPER  AND  LOWER  CUTOFF  FREQUENCIES  (TWO  ELEMENTS)  (RAD) 
CIS]  A 

CIS]  A  Y. KATZIH. I. A.F. .NOVEMBER  1987 


TULC1]«-0.4*PI 
TULC2  3*-0.2xPI 
n«-21 


hBS«-n  BSCOEFF  TUL 


7.799537305E*18  0.01284809274  *0.01443410434  *0.06997530689  *0.08163809137 
2.341241641E"17  0.1224571371  0.1632761827  0.05781641735  *0.113632834 
0.8  *0.1156328347  0.05781641735  0.1632761827  0.1224571371 
2. 341241641E*17  *0.08163809137  *0.06997550689  *0.01445410434 
0.01284809274  7.799S37305E*18 


▼HMBUNGCtm 
COJ  W*HAMMING  M;M;OIO 
Cl]  A 

C2]  A  GENERATING  THE  'N'  wide  HAMMING  WINDOW 

C41  010*0 

CS]  W*NpO 

t«3  •»t(H+2mf  CN+2]»/QD0 

C7]  A  EVEN  NUMBER  0?  SAMPLES  IN  WINDOW 

C8]  M*CN+2) 

C91  W*0.54*O.46xCOSCCCPIx23xul]+tn 
CIO]  W*(«W) ,W 
C11J  MS 

C12]  0D0:H*CS*l)+2 

C13]  W*0.34*0.46*C0SCCCPIx2Jxul)+N) 

C14J  W*(«W) ,liW 
CIS]  -*0 
CIS]  A 


C17]  A  THE  WINDOW  STARTS  AT  t»0  AND  HAS 
CIS  I  A  N<- NUMBER  OF  SAMPLES  CVECTOR] 

C 19  ]  A 


’M'  SAMPLES. 


C211  "  .NOVEMBER  1987 


N*U 

W*HAMMING  M 


vhanntngcoiv 
to]  W4-HANNING  M:M:OIO 
Cl]  8 

C2]  It  GENERATING  THE  'N'  WIDE  HANNING  WINDOW 

C  3  ]  It 

C4]  DIO*-0 

CS]  Wt-NpO 

CS]  -»CCN+2)*CrCN+2]])/ODD 

C7]  It  EVEN  NUMBER  OF  SAMPLES  IN  WINDOW 

CS]  H*-CN+2) 

C9J  W*-0.5*0.5*COSaCPIx2)*t.M]+N) 

CIO]  W4-C9W)  ,W 
Cll]  -*0 

C12  ]  ODD:M*-CN*l]+2 

C13  ]  W4-0.5*0.5i<C0SUCPI*2)*U!)+N) 

C 14 ]  WfCeWJ.llW 
CIS]  -*0 
CIS]  it 

C17]  It  THE  WINDOW  STARTS  AT  t«0  AND  HAS  'N'  SAMPLES. 
CIS]  It  N«- NUMBER  OF  SAMPLES  C VECTOR] 

CIS]  it 

C20]  a  E.KATZIR.I.A.F. .NOVEMBER  1987 
C21]  it 


N«-31 


W*- HANNING  N 


W 

2. 565338304E~3  0.0229303718  0. 0S282S69193  0.1206209387  0.1939470087 

0.2798029242  0.3746737339  0.4746754156  0.5757138888  0.6736526264 
0.7644820052  0.3444834595  0.9103817206  0.9594739058  0.9897649706  1 
0.9897649706  0.9594789058  0.9103817206  0.8444834595  0.7644820052 
0.6736526264  0.5757138888  0.4746754156  0.3746737339  0.2798029242 
0.1939470087  0.1206209387  0.06282669193  0.0229303718  2. 565338304E"3 


81 


VBARTLETTC03V 

toi  w* Bartlett  N:M:aio 
CU  A 

C21  A  GENERATING  THE  *H'  WIDE  BARTLETT  (TRIANGULAR)  WINDOW 
C31  A 

C4)  aio*-a 
C51  4«-N,  0 

C63  *( CN+2) T  (N+2 ) ) ) /ODD 

C7)  A  EVEN  NUMBER  OP  SAMPLES  IN  WINDOW 

C83  M«-CN+2) 

C9  )  W*-(A(l-(  IM)+M)  )  , l-( VM3+M 

CIO]  *0 

Cll]  ODDsM*CN+l)+2 
C121  W*-l-ClM)+K 
C13 )  R*C*W).UW 
C14]  -*0 

CIS]  A 

C16 ]  A  THE  WINDOW  STARTS  AT  t*0  AND  HAS  'N*  SAMPLES. 

C171  A  N«- NUMBER  OP  SAMPLES  (VECTOR] 

CIS]  A 

C191  A  E.KATZIR. I. A.P.  .NOVEMBER  1987 
C20]  A 


N«-31 


W*  BARTLETT  N 


W 

0.0625  0.123  0.1875  0.25  0.3125  0.375  0.4375  0.5  0.5625  0.625  0.6875  0.75 

0.8125  0.875  0.9375  1  0.9375  0.875  0.8125  0.73  0.6875  0.625  0.5625  Tt 
0.4375  0.375  0.3125  0.25  0.1875  0.125  0.0625 
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4-»  U»  U>  a.H  I 


VFREQRESCQ3V 
HW+M  FREQRES  h;aiO;N:H 
A 

A  SENERATING  FREQUENCY  RESPONSE  MAGNITUDE  AND  PHASE  FROM  TIME  DOMAIN 
A  COEFFICIENTS 

A 

aio«-o 

N+CM-ph)pO 
h+h.N 
H+PFT  h 

HW«-CM+2)TCXMAGN  H) 

PW«-(M+2  3  T  CXPHAS  H) 

-»0 

A  THE  FUNCTIONS  FFT ,  XMAGN . XPHASE  ARE  USED. 

A  IN  ORDER  TO  RECEIVE  A  BETTER  RESOLUTION  THE  COEFFICIENTS  ARE  2ERO- 
A  PADOET  TO  ’M’  SAMPLES. 

A  M*- NUMBER  OP  SAMPLES  INCLUDING  2ERO  PADDING  CUSE  ONLY  RADIX  TWO  NO. 
A  h*THE  FILTER  COEFFICIENTS  C VECTOR 3 
A  HW:  FREQUENCY  RESPONSE  MAGNITUDE  CM+2  SAMPLES) 

A  PW:  FREQUENCY  RESPONSE  PHASE  CM+2  SAMPLES)  CRAD.  ] 


Y.KATZIR.  I.  A.  F. NOVEMBER  1987 


hLP+21  LPCOEFF  C0.4XPI) 


HLP+128  FREQRES  hLP 


HLP 

0.9567807992  0.9619585014  0.9763182673  0.9965906571  1.018111173  1.03582952 
1.045421125  1.044271236  1.032124069  1.01125183  0.986091438  0.9624009 
0.9460845406  0.9419036347  0.9523233652  0.9767236974  1.011145443 
1.048645824  1.080227975  1.096202394  1.087755987  1.048461283  0.975462 
0.3701304604  0.7380554233  0.5883839646  0.4325882327  0.282876649 
0.1505027492  0.04424706035  0.03068671843  0.0731944037  0.08614789716 
0.07567134118  0.04996228114  0.01787457278  0.01248921845  0.0349417940 
0.04595168274  0.04486445606  0.03360579506  0.01595981635  3 . 4133U717E*3 
D. 02001949514  0.03043376941  0.03290460379  0.0275668884  0.01625216492 
1.963679943E~3  0.01185437487  0.02209076454  0.02662206686  0.024719091 
0.01713279563  5. 857126075E"3  6. 366569896E"3  0.01668743018  0.02278343 
0.02336291733  0.01342504363  9 . 228990525E"3  2. 014864013E"3  0.01265918 
0.02022375085 


■v'Cv.'C V't'C’t'  /  vl\  V"! 
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71DEAL7CD] V 

CO]  H^TYPE  IDEAL?  TCjQIO 
Cl]  A 

C2  3  A  GENERATING  IDEAL  FREQUENCY  RESPONSE 

C3  ]  A 

C41  0104-0 

CS1  H4-TCC0]<»0 

C6]  -*<Cx/CTYPE»’LP'33vCx/CTYPE»'HP'  )))/LP 

C7]  -*(Cx/(TYPE*'BP'  3]vCx/(TYPE»'BS'  )))/BP 

C 8 3  'IMPROPER  ARGUMENT. USE  ONLY:  LP.  HP  .  BP  .OR  BS  . 

C9]  -»0 

CIO]  LP:HC l( l  C  CTCCO  3+2  3  xTCC 13  3  *13 ]*T 

Cll]  H+CTCC01+23TH 

C12]  -Kk/CTYPE-'LP'  n/ continue 

C13]  H+-H 

C14  3  +CONTINUE 

C13]  BP:HCUClCCTCCO]+23x(TCCl]-TCC21  3  3  3  3  +  ClCCCTCCO]+23xTCC2]3  +  13  3]«-L 

C16]  H4-CTCC01+23TH 

C17]  ■•<  X/  CTYPE»  '  BP*  3  3  /CONTINUE 

C 18 1  h*— H 

C19  3  CONTINUE: H+H.OH 
C201  -*0 

C213  A 

C22]  A  'TCPE'+THE  TYPE  0?  FILTER  TO  BE  GENERATED  CLP. HP. BP, BS] 

C231  A  TCC 0 3 +NUMBER  OP  SAMPLES  CUSE  ONLY  RADIX  TWO  NO.  3 

C24  ]  A  TCC  13 ♦THE  CUTOFT  DIGITAL  FREQUENCY  FOR  LP  OR  HP  FILTER  OR 

C  25  3  A  THE  UPPER  DIGITAL  FREQUENCY  FOR  BP  OR  BS  FILTERS 

C26]  A  TCC 2  ]»THE  LOWER  DIGITAL  FREQUENCY  FOR  BP  OR  BS  FILTERS 

C27  ]  A  THE  FREQUENCIES  SHOULD  BE  A  FRACTION  OF  PI  CRAD.  ] . 

C28]  A 

C293  A  Y.KATZIR  .1. A.  F.  .NOVEMBER  1987 
C30]  A 


VFC0EFFCD3V 

COJ  h*-tl  FCOEFF  HW:QIO:H:a;SI2E 
Cl]  A 

C 2 3  A  CALCULATING  THE  'n'  FILTER  COEFFICIENTS  USING  IFTT 

C3  ]  A 

C4i  aio*-o 

C51  n«-n-l 

C6]  SI2E«-Un+21 

C7]  b«-SI2E*1 

C8J  H*-IFFT  HW 

C91  h«-Cl.a)tH 

CIO!  h«-mph 

C113  h«-CSI2ETC<Mi]).h 

C12]  -*0 

C131  A 

C 14 1  A  THE  COEFFICIENTS  ARE  FOR  CAUSAL  NONRE CURSIVE  FILTER. 

CIS]  A  n«-OOD  NUMBER  OF  REQUIRED  COEFFICIENTS  (RECTANGULAR  WINDOW) 
CIS]  A  HW*-  THE  IDEAL  FILTER  SAMPLES  (RADIX  TWO) 

C17]  A 

CIS]  A  Y, KAT2IR. I. A.F. .NOVEMBER  1987 
C19  ]  A 


H 

0000000011111111111090000000000000000 

oooooooooooooooooooooooooooooooooo 

0000000000000000000000000000000011 

11111100000000 

n«-21 

h«-n  FCOEFF  H 


0.02511120991  0.04315473177  0. 03628637547  ”5.2900 75645E”3  ”0.07065358471 

”0.1217862264  ”0.1233903877  ”0.06216894475  0. 0397303U337  0.133830098 
0.171875  0.1338300989  0.03973030337  ”0.06216894475  ”0.1233903877 
”0.1217862264  *0.07065358471  ”6. 29007564SE”3  0.03628637547  0.0431547 
0.02511120991 
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UNCLASSIFIED 


PC  SOFTNARE  FO*  THE  TEACHING  OF  DIGITAL  SIGNAL 
PROCESSINGS)  NAVAL  POSTGRADUATE  SCHOOL  HONTEREV  CA 
V  KATZ I R  HAR  M 
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9SACFtD37 

CO]  R«-LAG  SACF  DAT ; DIO; MEAN :N 

C23  fl  SAMPLE  AUTOCORRELATION  FUNCTION  (UNBIASED) 
C  3 1  fl 

C  4 :  QIO«-K«-0 

(5)  N«-pDAT 
C  6 1  R«-LAGpO 

(7)  DAT<-DAT .  R 

C  3  3  LOOP:RCK3«-C  +  /0ATx(-K)«DAT)+N 
[9)  -»(LAG>K«*K+l)/LOOP 
CIO)  +0 

C12)  fl  Y.KATZIR.I.A.F. .JULY  1987 

C14)  r  DAT«-DATA  VECTOR  CONE  DIMENSIONAL) 

(153  fl  LAG«-NO.  OF  LAGS  (MAX.  300) 

C161  R 


Rxx<-100  SACF  GAUSD 


Rxx 

1.1816  "0.011604  "0.066841  "0.14886  "0.11086  "0.0S0135  "0.28244 

0.12983  0.011973  "2.0869E"3  "0.052813  "0.066688  0.1414  0.16079 
0.021921  "0.050052  "0.03314  "0.12319  "0.11111  "0.087361 
7.9253E"3  0.047914  0.082207  0.022012  "0.072971  0.047907 
"0.028347  0.03337  "0.091018  0.11214  4.3453E"3  "0-020604  0.08519 
"0.13963  0.038409  "0.059126  "0.023359  "0.04759  0.03288  0.23321 
"0.021454  0.051991  "0.049765  *0.025113  0.059572  "0.11704  0.11551 
"0.11528  "0.058985  "0.038291  "0.077532  0.11274  0.10221  0.18426 
"0.09541  "0.14171  8. 3982E~3  "0.10611  "0.13464  0.037944  0.035897 
0.13416  "0.081343  "0.028463  0.13128  0.026911  0.011188  "0.048038 
0.048317  "0.023835  "0.14964  "0.071093  0.084411  0.04223  0.055516 
0.08221  "0.011808  "0.03549  "4.7687E"3  "0.011637  0.010587 
>.3499E"3  0.033257  "0.065615  "0.036728  0.038447  0.014501 
0.029666  0.042598  "0.096824  5.4581E"3  0.052421  0-047676 
"0.044033  "0.049402  0.016731  "1.8189E"3  6.5336E"3  "0.028012 
0.015022 
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tmeaned:v 

CO:  M«-MEAN  F 
£13  A  CALCULATE  THE  MEAN 
C23  M*-(  +  /*/F)+*/pF 

C33  A 

[43  A  SOURCE: MDSP  W. SPACE 
C53  A  F*-DATA  VECTOR 
[63  A  M:THE  MEAN  CNUMBER3 
[73  A 


X«-1100 


X 


3  4  5  6 

7  ( 

1  9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

XX*- MEAN  X 
XX 

50.5 

x*-x-xx 

X 


”48.5 

*47.5  "46.5 

'45.5  "44.5  "43.5 

"42.5 

"41.5 

”40.5 

"39.5 

"38.5 

“37.5  "36.5 

"35.5  "34.5  *33.5 

"32.5 

"31.5 

"30.5 

"29.5 

”28.5 

”27.5  ”26.5 

"25.5  "24.5  "23.5 

"22.5 

"21.5 

"20.5 

"19.5 

*18.5 

*17.5  "16.5 

"15.5  "14.5  "13.5 

"12.5 

“11.5 

“10.5 

”9.5  "8.5 

~7.5  "6.5  '5.5  ”4. 

5  “3.5  "2.5  "1.5  ' 

'0.5  0. 

5  1.5 

2.5  3. 

5  4.5 

5.5 

6.5  7.5  8.5  9.5  10.5  11.5  12.5  13.5  14.5  15.5  16.5  17.5  18.5 

19.5  20.5  21.5  22.5  23.5  24.5  25.5  26.5  27.5  28.5  29.5  30.5  31.5 

32.5  33.5  34.5  35.5  36.5  37.5  38.5  39.5  40.5  41.5  42.5  43.5  44.5 

45.5  46.5  47.5  48.5  49.5 


XX*- KEAN  X 


Jj 


■iju*. 


vcovcaiv 

x«-cov  p.;nro;N:i 

A 

A  GENERATE  TOEPLITZ  COVARIANCE  MATRIX  FROM  COVARIANCE  FUNCTION 

A 

aio«-i«*o 

N«-pR 

X«-CN.N)p0 

R«-(*liR).R 

LOOP: XI I:  3*-(-N)fC-I)iR 
-»(N>I«*I*1)  /LOOP 

A 

A  SOURCE:  MDSP  W. SPACE 
A  R<-THE  CORRELATION  FUNCTION  (VECTOR) 


RCLl<-3  SACF  CLI 


RCL1 

7.3328  6.8337  6.4035 


CMCL1«-C0V  RCLI 


GMCL1 

7.3328  6.8337  6.4035 
6.8337  7.3328  6.8337 
6.4035  6.8337  7.3328 


\X7 


y  jf  \» 
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VLPFPCOJV 

CO]  a«-P  LPFP  DATA ; DIO ; R ; K : MDATA : M 
CIJ  A 

C 2]  A  GENERATE  THE  Pth  ORDER  LINEAR  PREDICTION  FILTER  PARAMETERS 

[  3  ]  A  AND  THE  PREDICTION  ERROR  VARIANCE 

C4]  A 

CS3  010«-0 

C63  a«-CP+l]pO 

C7]  M«-MEAN  DATA 

C3:  MDATA«-DATA-M 

C9]  R«-CP+1)SACF  MDATA 

CIO]  K«-COV  R 

Cll]  a«-l.a«— C11R28CP.P]  TK 
C12  ]  PERR«-+/CaxR) 

C 13  ]  -*0 

C 14  ]  A 

CIS]  A  THE  FUNCTION  SOLVES  NORMAL  EQUATIONS  USING  THE  SACF.COV.AND  MEAN 
C16]  A  FUNCTIONS. 

C17]  A 

C 18 ]  A  Y.KATZIR.I.A.F. .AUGUST  87. 

C 19  ]  A 

C20]  A  a:  FILTER  COEFFICIENTS 

C21]  A  PERR:  THE  ESTIMATE  PREDICTION  ERROR  VARIANCE 
C  22  ]  A  Pt-THE  FILTER  ORDER 
C23  ]  A  DATA«-THE  DATA  VECTOR 
C24 ]  A 


a«-3  LPFP  CLl 


a 

1  ”0.897346235  "0.017308669  "0.02113346634 


PERR 

0.9625826774 


flHjUWWjj  y?  m  w  mro*  u  u  h  uv  wi[u  .  w  .*  w  a  jv  ^ 

i 
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i 
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VCCCD]V 

CO]  Y«-X  CC  H:QIO;N;I :NN:SNN 
Cl]  A 

C2 ]  A  1-D  CIRCULAR  CONVOLUTION 
C  3  ]  A 

C4 1  qio*-i«-o 

C5]  SNN«-pCNN«-pH) 

C6]  N*-pX 
C7]  •♦CNXNN)  /ERROR 
C8]  -KSNNj<1) /ERROR 
C9]  Y«-NpO 

CIO]  LOOP:YCI]**/XxC-l+I)®«H 
Cll]  -KN>I«-I+l)/LOOP 
C12]  *0 

C13J  ERROR- 'ARGUMENTS  ARE  NOT  THE  SAME  SIZE  OR  1-0. 
C14 ]  -»0 

C15]  A 

C16]  A  Y.KATZIR. I. A.F. AUGUST  1987 
C17]  A 

C18]  A  ARGUMENTS  NEED  TO  BE  THE  SAME  SIZE  AND  1-D. 
C19 ]  A  X  AND  H«- ARGUMENTS  TO  BE  CONVOLVED  C VECTOR) 
C20)  A 


X 

3  2  1 

H 

3  2  1 


Y«-X  CC  H 


A 


Y 

13  13  10 

LONG 

1. 23  4  56789  10 


i 


X'CC  LONG 

ARGUMENTS  ARE  NOT  THE  SAME  SIZE  OR  1-D. 


A 


•N 


a 


s 

*< 


93 


y 
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VFCCCO]V 

Co:  Y«-X  FCC  H;QIO:N:NN;SNN;A 

C13  A 

C23  A  FAST  1-D  CIRCULAR  CONVOLUTION 

C3]  A 

C4  ]  QXQ<-0 

C5]  SNN*-pCNN«-pH] 

C6]  N«-pX 

C7]  -»CNN>805/STATM 
C8  ]  -KN/NN) /ERROR 
C9]  -»CSNNj‘1] /ERROR 
CIO]  A«-CN.N]p<#H 
Cll]  A«"(  -lN)®A 
C12]  X«-l»CA+.xXD 
C 13  ]  -*0 

C14]  ERROR: ’ARGUMENT  ARE  NOT  THE  SAME  SIZE  OR  1-D 
C15]  -*0 

C16]  STATM: ' PLEASE  USE  CC  OR  LCV  FUNCTIONS. ’ 

C 17  ]  -*0 

CIS]  A 

C19]  A  X. KAT2IR.T. A. F. .AUGUST  1987 
C20]  A 

C21]  A  ARGUMENTS  NEED  TO  BE  THE  SAME  SIZE  AND  1-D. 
C221  A  LENGTH  OF  ARGUMENT  VECTOR  NEED  TO  BE  SHORT. 
C  23  ]  AX  AND  H«- ARGUMENT  TO  BE  CONVOLVED 
C24]  A 


3  2  1 


3  2  1 


Z«-X  FCC  H 


13  13  ro 
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f 
r 

c2 


a 

4 


i 

y? 

.5 

> 

.  •* 

wm 

!f 

V 

ft 


4 

y. 

V, 

a 


VLCVCCI3V 

CO]  Y«-xi  Lev  Hl;OlO:Ni:M2:N:X2:H2:X;H 

Cl]  A 

C2]  A  1-D  LINEAR  CONVOLUTION 

C  3  ]  A 

C  4  1  0I0«-0 

C5]  Nl«-pXl 

CS1  N2«-pHl 

C7]  N«-N1+N2-1 

CS  ]  X2«-NTX1 

C9]  H2<-NTH1 

CIO]  Y+-X2  CC  H2 

Cll]  -*0 

C 12  ]  A 

C13 ]  A  LINEAR  CONVOLUTION  USING  CC  FUNCTION 
C 14  ]  A 

CIS]  A  Y. KAT2IR. I- A.F. .AUGUST  1987 
C16]  A 

C17]  A  X  AND  H«-THE  ARGUMENTS  TO  BE  CONVOLVED  C VECTORS 5 


3  2  1 


3  2  1 


Y+-X  LCV  H 


9  12  10  4  1 


LONG 

1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29 

30-31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  5253 

54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77 

78  79  80  81  82  33  84  85  86  87  83  39  90 

Y*-X  LCV  LONG 

Y 

3  3  14  20  26  32  38  44  50  56  62  68  74  80  86  92  98  104  110  116  122  128  134  140 
146  152  158  164  170  176  182  188  194  200  206  212  218  224  230  236  242  248 

254  260  266  272  278  284  290  296  302  308  314  320  326  332  338  344  350 3C4 

362  363  374  380  386  392  398  404  410  416  422  423  434  440  446  452  458  464 

470  476  482  488  494  500  506  512  518  524  530  536  269  90 


VFLCVCOIV 

COI  Y«-X1  FLCV  H1:0I0;N1;N2 :N;XZ ;HZ;X:H 
CXI  A 

C2I  A  PAST  1-D  LINEAR  CONVOLUTION 

C3 1  A 

C4  3  010«-0 

C5I  Nl«-pXl 

C6I  N2«-pHl 

C7I  N«-N1+N2-1 

f  81  XZ*-NTX1 

C9I  HZ«-NTH1 

CIO 3  Y«-XZ  FCC  HZ 

C113  -*0 

CX2I  A 

C133  A  FAST  LINEAR  CONVOLUTION  USING  FCC  FUNCTION 
C14  3  A 

CIS!  A  Y. KATZIR. I. A.F. .AUGUST  1987 
C16I  A 

C173  A  X  AND  H«-THE  ARGUMENTS  TO  BE  CONVOLVED  C VECTORS I 
C18I  A 


3  2  1 


3  2  X 


Y*-X  FLCV  H 


Y 

9  12  10 


Y«-X  FLCV  LONG 

PLEASE  USE  CC  OR  LCV  FUNCTIONS. 
VALUE  ERROR 
FLCVC103  Y«-XZ  FCC  HZ 
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f: 

K 
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t 


VFFTCCJ3V 

C03  z«-FFT  x;k;n:r:aiO 
Cl]  dio+o 

C  2  ]  -»Cr«-(2€ppX)A2«lTpX]f8 

C3  ]  •*Cl«pk]T7-r<-25>k«'l0.5xpx 

C4]  'NOT  A  COMPLEX  VECTOR' 

C5]  -»0 

C  6]  -»8.pX«-X.C  "0.530 

C7]  x«-*Ck.23px+2 

C83  -»C20ik«-  "1  2  4  8  16  32  64  128  256  512  1024  2048  4096  8192  16384  3276 

C 9 ]  -»Cllpz«-X]l  0  0 

CIO]  'NOT  A  POWER-OF-2  ORDER' 

Cll]  -+-UEX  'Z' 

Cl'.  ]  Z«-C1  2  ,kp23px 

C 13  ]  -»ki  0  28  23 

C14  ]  -»<2<aNC  'wt'Uie 

CIS]  -*(k«ppx«-wt]T20 

C16]  n«-o,-«14.x«-2oCo2+n)xtn+4 

C 17 ]  vt«-x«-Ckp23p(x,n3  .n.-x 

C 18  ]  -*20 

C19]  x«-«ckp’xC: ’0]' 

C20]  n«--/C23z 

C  213  k*pz«-+/C23z 

C22]  •*(4<k«-ppz«-z.C0]C-/cl]nxkpx]  .C0.5]*/Cl]nxkpex]ti9 

C  23 ]  n«--/C2  3z 

C24  ]  z«-+/C2]z 

C25]  nC:0:13«--nC;0:i3 

C26]  nC;  1  0  ;13«-nC :  :13 

C27]  z«-z,C0]n 

C 28 ]  z«-«*C*/z]  ,tn«-0]-/z 

C29]  -Ur 

C30]  n«-”i®«z 

C 3 1  ]  nCl;  ]«~nCl;  ] 

C32  3  x«-n-z 

C  3  3  ]  n«-"ltpz«-z+n 

C  34  ]  k<-0,®Hr«-2o(o+n)xvn+2 

C35]  r«-.Cpz]pk.r,-r 

C36]  x^-c+zxxr]  ,C"0.5]-/rxex 

C37 ]  z*Ci*x) ,z-x 

C  38  3  -*0 

C393  ft 

C40]  ft  APL- STYLE  FAST  FOURIER  TRANSFORM 

C41J  ft  x«-CAN  BE  EITHER  REAL  CA  VECTOR]  OR  COMPLEX  C2-ROW  MATRIX] 

C 42 3  ft  OPTIMIZED  FOR  SPEED 

C 43 ]  ft  ZsTHE  OUTPUT  IS  COMPLEX  CREAL  AND  IMAGE  PARTS]. 

C 44 3  ft  SOURCE : PAUL  PENFIELD . JR. 6/ 13 /79 
C 45 3  A 


V 


v 


*s 
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*  ^ 


viFFTCtm 

CO]  z*-IFFT  x:k;n;r;OIO 

ci]  aio«-o 

C  2  ]  -*Cr<-(2«ppX)A2«lTpx)T8 

[3  ]  -*Cl<pk)T7-r«-25>k«-l.0.5xpx 

C 4 ]  'NOT  A  COMPLEX  VECTOR' 

C  S  ]  -*0 

C6]  -*8  ,  px«-x.  C  “0. 5] 0 

C7]  X«-«Ck.2)pX+4 

[8]  -»(20ek«-  *1  2  4  8  16  32  64  128  256  512  1024  2048  4096  8192  16384  3276 

C 9  ]  -*(14pz«"X]  i  0  0 

CIO]  'NOT  A  POWER-OF-2  ORDER' 

till  -*-QEX  '2' 

C 12 ]  z«-(l  2  ,kp2)px+n 

C13  ]  -*ki  0  28  23 

C 14  ]  -K240NC  'Wt'U16 

C 15 1  -»Ck«ppX«-wt)  T20 

C16]  n«-0.-<»lix«-2o(o2+n)xi.n+4 

C 17]  wt»-x«-(kp2)pCx.n)  .n.-x 

1 18  ]  -*20 

ci9 ]  x«-«(kp'xc ; : '0] ' 

[20]  n«--/C2]2 

C21]  k«-pZ*-*/C2]Z 

C22 ]  -»(4<k«-ppz«-Z.C0]C+/Cl]nxkpx]  .  CO.  5W  C l]nxkpe-x)T19 

C23 ]  n*~/C2]z 

C24  ]  Z*-+/C2]2 

[25]  nC :l:l]«—nC :1:1J 

C26]  nC;  1  0  :l]*-nC :  ;1] 

[27]  z«-z.CO]n 

C28  ]  z«-«C*/z)  .Cn«-0]-/z 

C29  ]  -M.r 

C30]  n«-”l®<6z 

C  31]  nCl; ]«— nCl; ] 

C32]  x«-z-n 

[33]  n«-*lTpz«-z+n 

[  34  ]  k«-Ot«Hr<-20(o+n]xi.n+2 

C351  r«-(pz)pk.r.-r 

[36]  x«-(-/xxr] ,  C  ”0.  SWxxer 

C37  ]  Z«-(2+X),2-X 

[  38  ]  -*0 

C 39 ]  ft 

C40]  n  APL-STYLE  INVERSE  FOURIER  TRANSFORM 

[ 41 ]  ft  x«-CAN  BE  EITHER  REAL  CA  VECTOR]  OR  COMPLEX  C2-ROW  MATRIX) 

[ 42 ]  ft  OPTIMIZED  FOR  SPEED 

[43]  ft  SOURCE: PAUL  PENFIELD, JR.  6/13/79 

[44]  ft 
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1  .«  - 


A 


1 


X 

1111111111111111 

X«*FFT  X 

y 

16  000000000000000 

oooooooooooooooo 


r, 

r. 

h1 

f 

* 


A 


y 

16  000000000000000 

oooooooooooooooo 

XX*-IFFT  y 

XX 

1111111111111111 

oooooooooooooooo 

xx«-xxco:  J 

XX 

1111111111111111 
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^  ^  ^  JT  -A  J  s- -  +m 


--  V-.-.  «T.  -r.-\  S- 


PTC 


U  H1  '-l1  m1  l'-,Mi  *.■  l re » <■  cn.wx^  v- 


VPSECO]* 

to:  INW*-ZP  PSE  DATA ; QXO ;N ;  XN ;  XW 
Cl]  A 

C2]  A  COMPUTING  THE  POWER  SPECTRUM  ESTIMATION  CPERIODOGRAM) 

C  3  ]  A 

C  43  010«-0 

C5]  N*-(pDATA)+ZP 
C  6 1  N«-2*rC2«N) 

C7]  XN^DATA, C  CN-CpDATA3 ]p0] 

C8]  XW«-FFT  XN 
C9J  INW«-C+i‘CXW*2))+N 
CIO]  -*0 
Cll]  A 

C 12  ]  A  COMPUTING  THE  POWER  SPECTRUM  ESTIMATION  C  PERI 0 DOGRAM 3  USING  FFT 
C13  ]  A  FUNCTION 

C 14 1  A 

CIS]  A  Y.KATZIR  .I.A.F. .AUGUST  1987 
C16]  A 

C17]  A  DATA*THE  DATA  VECTOR  Cl-D) 

C18]  A  ZP«- LENGTH  OF  ZERO  PADDING  FOR  SMOOTH  PLOT 
C19]  A  INW.-THE  OUTPUT-SXCu] 

C 20 ]  A 


VPSEB  C □ ] 7 

CO]  INBWoSEG  PSEB  DATA:DIO:M:MM;XW;XN:L:X 
Cl]  A 

C21  A  COMPUTING  THE  POWER  SPECTRUM  ESTIMATION  USING  BARTLETT  WINDOW 

C3  ]  a 

C4  ]  010«-0 

CS]  M«-l  (  CpDATA]+SEG) 

C6]  XN«-CM,SEG]  pO 
C7]  ZP*-256-SEG 

C3  ]  MM«-CHpXN)*ZP 
C9]  XW«-CM.MM]pO 
C10J  L«-K«-0 

Clll  LOOP1:XNCKuSEG]*-DATACCSEGxK]*t.SEG] 

C 12  ]  •*(M>K«-K*1]  /LOOP1 

C13  ]  LOOP2:XWCL;]«-ZP  PSE  XNCL;  ] 

C 14  ]  •*CM>L*-L*1)  /LOOP2 

C15  ]  INBWM  +j‘XW)*M 
C 16]  -»0 

C 17  ]  A 

C18 ]  A  Y.KATZIR.I.A.F. .AUGUST  1987 
C19  ]  A 

C20]  A  SEG«-LENGTH  OF  SEGMENTS 
C21]  A  DATA*- THE  DATA  VECTOR  CI-D) 

C22 ]  A  USING  THE  PSE  FUNCTION  C *FFT) . 

C23 ]  A  THE  FUNCTION  PADS  THE  SEGMENT  WITH  ZEROS  UP  TO  256  POINTS. 

C24 ]  A 


»^r 


vxsqrtcu]? 

CO]  X2*-XSQRT  XR 
Cl]  A 

C2]  A  COMPLEX  SQUARE  ROOT  CXR*0.5) 

C  3  ]  A 

C4]  XZ«-C0.5-0a.>XR]*(2xCC+j‘XRxXR]*0.5)*I 

C5]  X2«*(0a.>  1  0  (‘XRJeXZ.CDIO-O.SKO  1  * 


X 

0.5  0.366 

XSQRT  X 

0.86602  0.49999 


VXTIMECajV 
CO]  XZ«-XL  XTIME  XR 
tl]  A 

1 2 ]  A  COMPLEX  TIMES  CXLxXR) 

1 3  ]  A 

C4]  XZ«-(C1  0  /XL)x  1  o  (‘XR)-C0  1  /XL)x  0 
C5]  X2«-XZ.CaiO](Cl  0  /XL]x  o  1  (*XR)  +  C  0  1 


X 

0.5  0.366 


1  0  ♦ . xXR) *0. 5 
.xXR)+2xXZ+X2=0 


1  (‘XR 

/‘XL)  x  1  0  (‘XR 


X  XTIME  X 
"0.49996  0\ 866 


▼XCONJCQ]? 

CO]  XZ«-XCONJ  XR 
Cl]  n 

C2]  n  COMPLEX  CONJUGATE  C+XR] 

C  3  ]  n 

C4]  XZ«-C1  0  /XRJ.COIO]-  0  1  /XR 


0.5  0.866 


XCONJ  X 
0.5  ~0. 866 


VXEXP0CQ]V 
CO]  XZ«-XEXP0  XR 

cu  « 

C2]  A  COMPLEX  EXPONENTIAL  C*XR] 

C  3  3  A 

C4]  XZ«-CCpXR]po  1  o  /XR] *  2  1  o.o  0  1  o.xXR 


0.5  0.866 


XEXPO  X 
1.0682  1.2559 


. .HLK.^.V k ■ 


<iV 


) LOAD  0  EC4440 

0  EC4440  SAVED  2/23/1988  14:26:16 


NAVAL  POSTGRADUATE  SCHOOL 


MONTEREY.  CALIFORNIA 


VERSION  1.0 


ECZ-4440  WORKSPACE 

THESIS  BY 
Y.KATZIR.I.A.F. 

ADVISOR: PROF.  C.W.  TKERRIEN.  . 

SEPTEMBER  1987 


)FNS 

BARTLETTC  BARTLETTR  CC  CC2D  CHEB  COEFF  COS  FCC 

FCC20  FFT  PFT2D  FLCV  HAMMINGC  HAMMINGR  HANNINGC 

HANNINGR  IFFT  IFFT2D  ILPFILX  INFORMATION  LCV  MCCLEL  MEAN 

PI  PROTFILT  RECTANR  SHIFT  SHIFT2D  SIN  TRANS FNC  UNSHIFT 

UNSHIFT20  XCONJ  XEXPO  XMAGN  XPHAS  XSQRT  XTIME 


Ww 


55W95! 


F»»  *1 


VCC2D(03V 

C03  ¥«-X  CC2D  H;I1;I2;DI0:N1;N2  :T1;T2 

C13  a 

C23  A  2-D  CIRCULAR  CONVOLOTION 
[33  a 

[43  -»CCA/(pX3«pH3A2«ppX3/BEGIN 
[53  0*' IMPROPER  ARGUMENTS' 

[63  -*0 

[7  3  BEGIN :  0IO«-0 
[8  3  Nl**ltpX 

[9  3  N2«-lipX 

[103  Y«-CN2,N13puO 
[113  Il«-0 

[123  L00P1:T1«-C-CI1*1)  3*<#H 
[13  3  I2«-0 

[143  LOOP2:T2«-(-CI2  +  l)  )««T1 
[153  Y[I2:I13«-+/*/XxT2 

[16  3  ■*CN2>I2«-I2*13  /LOOP2 

[173  -*<N1>I1*-I1+13  /LOOP1 

[183  -*0 

[193  A 

(203  A  THE  INPUTS  NEED  TO  BE  2-D  AND  THE  SAME  SIZE. 
[213  A 

[223  A  T.KATZIR  I. A. F. .SEPTEMBER  1987 
[233  A 
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VFFT2DCD3V 

t03  Y*ITT2D  X:OIO:NROW;NCOL:I:J;G 
C13  A 

C23  A  2-DIMENSIONAL  ITT 

C3 :  a 

C4  3  £3IO«-0 

C53  A  CONVERT  REAL  ARRAY  TO  COMPLEX  ARRAY 
C63  I«-pCpX) 

C71  -*(1*3 3  /CONTIN 

C83  NROW«-pXC:13 

C9  3  NCOL«-pXCl:J 

CIO  3  G«-C2.NROW.NCOL3pO 

CU3  GC0::3«-X 

C123  X«-G 

C13  3  A  COMPUTATION  OF  FFT 

C143  CONTIN: NROW^-pXCl;;  13 

C153  NCOL«-pXCl;l:  3 

CIS  3  G*-XxO 

C17  3  Y«-G 

C18  3  IKT«-0 

C193  COLLOOP:GC  ;  :J3«-FFT  XC  :  ;J3 
C203  -»(NCOL>J«-J*13/COLLOOP 
C213  ROWLOOP:YC:I:3«-FFT  GC;I;3 
C223  -*CNR0W>I«-I*13/R0WIX30P 

C233  A  IF  OUTPUT  ARRAY  IS  REAL  VALUE, 2ND  PLANE  IS  DELETED. 
C243  GCO;  ;3«-Yto:;3*YCl:;3 
C253  GC1:  ;3«-YCl;  :3 
C263  !*♦/♦/♦/ CGfiY) 

C273  -*(1*03  /STPNOW 

C283  -*0  ' 

C293  STPNOMiYt-Cl.NROW.NCODTY 
C303  Y«'CNROW,NCOL3pY 
C313  A 

C323  A  FFT  FUNCTION  HAVE  BEEN  USED. 

C33 3  A 

C343  A  Y. KATZ IR. I. A. F.  .SEPTEMBER  1987 
1 3  5  3  A 
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VIFFT2DC01V 

COl  Y«-IFFT2D  X:aiO:NROW;NCOL;I:J;G 
Cll  A 

C21  A  2 -DIMENSIONAL  INVERSE  EFT 
C31  A 

C4 1  aia«-o 

C51  A  CONVERTING  REAL  ARRAY  TO  COMPLEX  ARRAY 

C61  I«-pCpXD 

C71  -*  ( I  *  3 )  /  CO  NT  IN 

C8  3  NROW«.pXC:ll 

C91  NCOL*-pXCl:I 

C101  G*-(2.NROW.NCOLlpO 

cm  gco:  :i<-x 

C123  X«-G 

C131  A  COMPUTATION  OF  XFFT 

C 14 1  CONTIN:  NROW«-pXCO ;  ;11 

C1S1  NCOL«-pXCO  ;1;  1 

C161  G*"X*0 

C171  Y«-G 

C181  W«-0 

C191  COLLOOP:GC;:Jl«-IFFT  XC  ;  ;J1 
C20 1  *(NCOL>J«-J+ll/COLLOOP 
C 211  ROWLOOP:YC  ;!;]«■  IFFT  GC;X:1 
C22 1  •*(NROW>I«-I+ll  /ROWLOOP 

C231  A  IF  OUTPUT  ARRAY  IS  REAL  VALUE, 2ND  PLANE  IS  DELETED. 

C241  GCO ; ; l^YCO; ;1+YC1;;1 

C251  GC1;;1«-YC1;;1 

C261  I«-+/  +  /-*'/CG*Yl 

C271  -»  C I  ■  0  )  /  STPNQW 

C281  -»0  . 

C291  STPNOW:Y«-Cl,NROW,NCOL)TY 
C301  Y^CNROW.NCOLlpY 
C311  A 

C321  A  IFFT  FUNCTION  HAS  SEEN  USED. 

C33 1  A 

C341  A  Y.EATZXR.I.A.F. .SEPTEMBER  1987 
C351  A 


00000000 

00000000 

Y«-FFT2D  HI 


H2«-IFFT2D  Y 
H2 

0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
oooooo 
0  0  0  0  0  0 
oooooo 
oooooo 

XI 


V«-FFT2D  XI 
Y 

[ 

I 


X2«-IFFT2D  Y 


z 

0000000000000 

0000000000000 

0000000000000 

0000000000000 

0000000000000 

0000000000000 

0000000000000 

0000000000000 

0000000000000 

0000000000000 

ooooooooooooo 

0000000000000 

ooooooooooooo 

ooooooooooooo 

ooooooooooooo 

ooooooooooooo 

Z1«-FFT2D  Z 
Z2«-XMAGN  Z1 

Z2 


8.5 

7.2 

5.3 

3 

0.7 

1.2 

2.5 

3  . 

2.5 

0.7 

3 

5.3 

7.2 

8.5 

8.1 

6.9 

5 

2.3 

0.67 

1.2 

2.4 

2.8 

2.4 

0.67 

2.8 

5 

6.9 

8.1 

6.9 

5.8 

4.3 

2.4 

0.57 

1 

2 

2.4 

2 

0.57 

2.4 

4.3 

5.8 

6.9 

5 

4.3 

3.1 

1.8 

0.41 

0.73 

1.5 

1.8 

1.5 

0.41 

1.8 

3.1 

4.3 

5 

2.8 

2.4 

1.8 

1 

0.23 

0.41 

0.85 

1 

0.85 

0.23 

1 

1.8 

2.4 

2.8 

0.67 

0.57 

0.41 

0.23 

0.055 

0.097 

0.2 

0.23 

0.2 

0.055 

0.23 

0.41 

0.57 

0.67 

1.2 

1 

0.73 

0.41 

0.097 

0.17 

0.35 

0.41 

0.35 

0.097 

0.41 

0.73 

1 

1.2 

2.4 

2 

1.5 

0.85 

0.2 

0.35 

0.72 

0.85 

0.72 

0.2 

0.85 

1.5 

2 

2.4 

2.8 

2.4 

1.8 

1 

0.23 

0.41 

0.85 

1 

0.85 

0.23 

1 

1.8 

2.4 

2.8 

2.4 

2 

1.5 

0.85 

0.2 

0.35 

0.72 

0.85 

0.72 

0.2  . 

0.85 

1.5 

2 

2.4 

1.2 

1 

0.73 

0.41 

0.097 

0. 17 

0.35 

0.41 

0.35 

0.097 

0.41 

0.73 

1 

1.2 

0.67 

0.57 

0.41 

0.23 

0.055 

0.097 

0.2 

0.23 

0.2 

0.055 

0.23 

0.41 

0.57 

0.67 

2.8 

2.4 

1.8 

1 

0.23 

0.41 

0.85 

1 

0.85 

0.23 

1 

1.8 

2.4 

2.3 

5 

4.3 

3.1 

1.8 

0.41 

0.73 

1.5 

1.8 

1.5 

0.41 

1.8 

3.1 

4.3 

5 

6.9 

5.8 

4.3 

2.4 

0.57 

1 

2 

2.4 

2 

0.57 

2.4 

4.3 

5.3 

6.9 

8.1 

6.9 

5 

2.3 

0.67 

1.2 

2.4 

2.8 

2.4 

0.67 

2.3 

5 

6.9 

8.1 

112 


VSHIFTC03V 
CO]  Y ♦■SHIFT  X 
ci:  a 

[2  3  A  SHIFTING  PERIODIC  DATA  TO  AN  INTERVAL  STARTING  AT  THE  ORIGIN 
[33  A 

C 4  3  Y*-C  L  (  “l+pX3-i-2  )®X 
[53  -»0 

[63  A 

C71  A  SOURCE :M-D  SIGNAL  PROCESSING  W.S. 

C81  A 


000010000 


X1<-SHIFT  X 


1000  0-0000 


VUNSHIFTCOIV 
C01  Y«- UNSHIFT  X 

Cll  A 

C21  A  SHIFTING  PERIODIC  DATA  TO  AN  INTERVAL  CENTERED  AT  THE  ORIGIN 
[33  A 

C  4 1  Y<-C-LC~l+pX3+2  3«X 

[53  *0 

C61  A 

C73  A  SOURCE: M-D  SIGNAL  PROCESSING  W.S. 

CS3  A 


100000000 

X2«-UNSHIFT  XI 


000010000 


v; 


KhK' 


T.V7.V 


vshift2dcu3? 

Y«-SHIFT2D  X 
ft 

ft  SHIFTING  2-D  PERIODIC  DATA  TO  AN  INTERVAL  STARTING  AT  THE  ORIGIN 

ft 

I«-(t  C  "l*ltlipX)+2)«a  C  ■l*ltpX3  +  2  )*X 
-♦0 
ft 

ft  SOURCE: M-D  SIGNAL  PROCESSING  W.S. 

ft 


VUNSHIFT2DC03V 
Y«-UNSHIFT2D  X 
ft 

ft  SHIFTING  2-D  PERIODIC  DATA  TO  AN  INTERVAL  CENTERED  AT  THE  ORIGIN 

A 

X<-C-U"I+lflipX)+2  3®C-U"I+lTpX)  +  2)eX 
-*0 
ft 

ft  SOURCE: M-D  SIGNAL  PROCESSING  W.S. 

ft 


i 


9 


0  0  0  0  0  0  0 

0  0  0  0  0  0  0 

0  0  0  0  0  0  0 

0  0  0  0  0  0  0 

0  0  1  0  0  0  0 

3  0  0  0  0  0  0 

3  0  0  0  0  0  0 

3  0  0  0  0  0  0 

3  0  0  0  0  0  0 

Z1«-SHIFT2D  2 


100000000 

000000000 

000000000 

000000000 

000000000 

000000000 

000000000 

000000000 

000000000 


Z2«-ONSHIFT2D  Z1 
22 

0  0  0  0  0  0 
0  0  0  0  0  0 
oooooo 
0  0  0  0  0  0 
0  1  0  0  0  0 
oooooo 
oooooo 
oooooo 
oooooo 


M90SHHB 


VILPFILTCO]* 

CO]  W«-R  ilpfilt  P:DI0;I:J;RR:RB:PB 
C13  a 

C23  A  IDEAL  LP  FILTER  C CIRCULAR  BASE]  FREQUENCY  RESPONSE 

C3  ]  A 

C4]  OIOO 

C5]  W«-CP.PDpO 

C  6]  I«-J*-0 

C7]  PB«-CP+2)-l 

C  8  ]  RB«-CPxR]  +  C2xPI) 

C9]  LOOP:ER«-CCCI-PBD*2)+CCJ-PB)*2D3*0.5 

CX03  ■* C RR>RB ]  / ZERO 

C113  WCI;J]«-1 

[12  3  ZERO :-»CP>J«-J+ ID /LOOP 

C13]  J«-0 

C14]  -»CP>I+-I*lD/LOOP 
CIS]  -*0 
[16]  A 

C17]  A  R<-THE  RADIUS  OF  THE  SASE  CGAIN*1D 
C18]  A  P*-THE  DIMENSION  OF  THE  MATRIX 
C19]  A  W:THE  IDEAL  FREQUENCY  RESPONSE 
C20]  A 

C21]  A  Y.KATZIR.I.A.F. ,  OCTOBER  1987 
C22]  A 


■  ii»  v  K.'*_<.txtv_w:'iLw  K-'.  x.”  x-  x-i  xi  x-  x-  v.  x-  x  •  xr. 


VHAMMINGRC 03V 

WR+-II  HAMMINGR  N:OIO;HW1 ;HW2 : J;I:ROW; COL:S 

A 

A  RECTANGULAR  BASE  HAMMING  WINDOW 

A 
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APPENDIX  D 

COMPUTER  ASSIGNMENTS  AND  SOLUTIONS 

This  appendix  provides  a  few  samples  of  computer  assignments  and  related 
solutions  using  the  software  package.  These  are  in  addition  to  the  examples  in 
Chapter  III.  The  computer  assignments  in  EC  3400  were  taken  from  FIRST  PRIN¬ 
CIPLES  OF  DISCRETE  SYSTEMS  AND  DIGITAL  SIGNAL  PROCESSING  by 
Professors  R.  D.  Strum  and  D.  E.  Kirk  [4].  The  computer  assignments  in  EC  3410 
and  EC  4440  were  taken  from  Professor  C.  W.  Therrien’s  homework  assignments. 


EC  3400  COMPUTER  ASSIGNMENT  [4] 

The  purpose  of  a  bandstop  filter  is  to  remove  the  effects  of  the  vibrations  in 
the  frequency  range  6  Hz  to  15  Hz  from  a  control  loop.  The  sampling  frequency  is 
fs=100  Hz. 

a.  Using  the  Fourier  Series  design  procedure,  calculate  the  filter  coefficients  with 
a  rectangular  window. 

b.  Plot  the  magnitude  portion  of  the  frequency  response  for  a  filter  having  thirty- 
one  coefficients. 

c.  Repeat  part  (b)  using  a  von-Hann  window. 


SOLUTION 

1.  By  using  the  function  DIGFREQ  the  digital  upper  and  lower  frequencies  are 
calculated  (Figure  D.l).  Figure  D.2  describes  the  ideal  frequency  response  to 
be  generated. 


TuL»-2flO 

TOLC13«-100DIGFREQ  15 
TULC13 
0.94243 

TULC13+PI 

0.3 

TOLC21«-100  0IG7RZQ  6 
T0T.C21 
0. 37699 

TULC23+PI 

0.12 


Figure  D.l.  Calculating  the  Digital  Frequency  Using  the  Function  DIGFREQ. 
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Figure  D.2.  The  Ideal  Frequency  Response. 

2.  By  using  the  function  BSCOEFF  31  causal  impulse  response  coefficients  axe 
calculated  (See  Figures  D.3  and  D.4). 

3.  The  function  HANNING  generates  31  von- Harm  window  samples  that  axe 
multiplied  by  the  filter  coefficients.  The  result  is  windowed  coefficients  (See 
Figures  D.5  and  D.6). 

4.  By  using  the  function  FREQRES  the  frequency  responses  with  and  without 
window  are  generated.  The  following  shows  the  APL  commands  which  gener¬ 
ate  the  frequency  responses,  and  Figures  D.7  and  D.8  show  the  plotted  results. 

HBS«— 128  FREQRES  hBS 
HBSW<— 128  FREQRES  hBSW 


hBS«-31  3SC0BT7  TOL 


*0.033694  *0.032361  *0.016485  *3.28428*4  *1.0218E“3  *0.01871  "0.037409 

*0.032854  7.33488*3  0.07206  0.12421  0.12619  0.063217  *0.042416  *0. 
0.32  *0.14034  *0.042416  0.063217  0.12619  0.12421  0.07206  7.35438*3 
*0.032854  *0.037409  *0.01371  *1.02188*3  *8.23428*4  *0.016483  *0.03 
*0.033694 


Figure  D.3.  Calculating  the  Impulse  Response  Coefficients 
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Figure  D.4.  The  Impulse  Response  Coefficients. 
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2.5653E”3  0.02293  0.062827  0.12062  0.19395  0.2798  0.37467  0.47468  0.5737 
0.67365  0.76448  0.34448  0.91033  0.95948  0.98976  1  0.98976  0.95943 
0.91038  0.84448  0.76448  0.57363  0.37571  0.47468  0.37467  0.2798  0.194 
0.12062  0.062327  0.02293  2.5653E"3 
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Figure  D.5.  Generating  the  Hanning  Window. 


EC  3400  COMPUTER  ASSIGNMENT  [4] 

A  txonrecursive  bandpass  filter  is  to  be  designed  using  the  IDFT  approach. 
The  ideal  characteristic  is  given  below: 


5ff  7T  ^  7T  57 r 

_?r  ~~6  ~2  °  2  T  ^ 

Use  255  frequency  samples  and  find  51  filter  coefficients  vising: 

a.  A  rectangular  window 

b.  A  Hamming  window 

c.  Plot  the  frequency  response  for  the  filters  of  parts  (a)  and  (b). 


SOLUTION 

1.  By  using  the  function  IDEALF  and  choosing  256  samples,  the  ideal  bandpass 
frequency  response  is  generated  (See  Figure  D.9). 

2.  By  using  the  function  FCOEFF  51  coefficients  are  calculated  (See  Figure 
D.10). 

3.  The  function  HAMMING  generates  51  Hamming  window  samples  that  are 
multiplied  by  the  filter  coefficients.  The  result  is  the  windowed  coefficients. 


V  AT.  X."WA.'JW--V-'-  . V.  . 
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TCCi:«-256 
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Figure  D.9.  Generating  the  Ideal  Frequency  Response  Using  255  Sam¬ 
ples. 
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"6.2133E"3  "4.642E“3  2.S82E*3  0.020547  "0.027498  4.703SE"3  8.7724E"3  4.8s*3 
"3.5166E"!  "0.023638  0.040397  *0.010498  "0.013214  “5.0631E"3  0.011 
0.036326  "0.069496  0.023239  0.024298  S.171SE"3  "0.029088  *0.077913 
0.21179  "0.12868  "0.16131  0.32813  "0.16131  "0.12868  0.21179  "0.077 
"0.029088  5. 1718E“3  0.024298  0.023239  "0.069496  0.036326  0.011138 
"3. 0631E“3  "0.013214  "0.010498  0.040397  “0.025638  "5.5166E"3  4.885E-J 
8. 7724E“3  4. 7036E"2  “0.027498  0.020547  2.582E“3  "4.642E"3  "6.2133E"3 


Figure  D.10.  Calculating  the  Impulse  Response  Coefficients. 


IM-HAMMING  31 


KB  P—128  FREQRES  h 
HBPW—128  FREQRES  hW 


Figure  D.ll.  Generating  the  Hamming  Window  and  Calculating  the 
Frequency  Response. 
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EC  3410  COMPUTER  ASSIGNMENT  [5] 


On  a  floppy  disk  you  will  find  4  data  sets  called 

SLO .DAT 
SL1.DAT 
SL2.DAT 
SL3.DAT 

Each  of  these  represents  a  random  signal  with  512  time  points.  The  format  of  each 
file  is  identical.  You  can  look  at  the  file  with  an  editor.  The  first  line  contains  the 
number  512  displaying  the  number  of  points  in  the  data  set.  On  the  remaining 
lines,  having  four  numbers  per  line,  you  will  find  the  floating  point  values  of  the 
signal. 

1.  Transfer  these  files  to  your  APL  WS  using  the  UTILITY  function  GETDATA. 

2.  Plot  eacft  of  the  signals  versus  time  and  submit  the  plots  with  your  assign¬ 
ment.  Tell  whether  you  think  each  of  these  signals  seems  to  be  more-or-less 
uncorrelated,  positively  correlated,  or  negatively  correlated. 

3.  Compute  the  mean  of  each  signal  by  averaging  in  time  and  write  down  your 
result  for  each  signal.  It  is  reasonable  based  on  your  plots? 

4.  Subtract  the  mean  from  each  signal  and  compute  and  plot  the  sample  corre¬ 
lation  function.  Attach  the  plots  to  your  assignment.  Was  your  guess  about 
the  correlation  of  these  signals  in  part  2  correct? 


SOLUTION 

1.  Using  the  function  GETDATA  and  the  following  statements,  the  data  sets 
SLO,  SL1,  and  SL3  are  transferred  to  the  APL  workspace: 

SLO*— GETDATA  ’ SLO .DAT’ 

SL1<— GETDATA  ’SL1.DAT 
SL2<— GETDATA  ’ SL2.DAT 
SL3*— GETDATA  ' SL3.DAT 


2.  By  using  the  X-Y  plotting  procedure  in  STATGRAPHICS  the  data  plots  are 
generated.  (See  Figures  D.15,  D.16,  D.17,  and  D.18.)  Another  kind  of  plot 
can  be  obtained  by  using  the  STATGRAPHICS  Time  Series  Procedure.  The 
results  are  shown  in  Figures  D.19,  D.20,  D.21,  and  D.22. 
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Figure  D.17.  SL2  Data  Set. 
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Figure  D.18.  SL3  Data  Set. 

3.  From  the  plotting  of  each  data  set,  it  can  be  seen  that: 

a.  SLO  changes  randomly  with  little  consistency,  this  is  characteristic  of  a 
white  noise  sequence. 

b.  SLl  varies  more  slowly  and  stays  stable  for  a  while  with  positive  values 
and  then  changes  more-or-Iess  smoothly  to  negative  values.  This  indicates 
that  it  is  a  positively  correlated  data  set. 

c.  SL2  is  very  nearly  like  a  white  noise  sequence. 


153 


3.  3 
2.  S 
1 . 5 
O.  3 
-O .  3 
-1. 3 
-2.  3 


O  ICO  200  300  400  300  600 

Tim*  C  N3 


Figure  D.19.  SLO  Data  Set  Using  Time  Series  Procedure. 


Figure  D.20.  SL1  Data  Set  Using  Time  Series  Procedure. 

d.  SL3  changes  at  almost  every  step  from  negative  to  positive.  This  is  char¬ 
acteristic  of  a  negatively  correlated  data  set. 

4.  Calculate  the  mean  of  each  data  set  using  the  MEAN  function.  Then  subtract 
the  mean  of  each  data  set  from  the  data  set  in  order  to  achieve  a  zero  mean  data 
set.  Using  the  following  instructions,  calculate  the  first  150  autocorrelation 
lags  for  each  data  set: 


Rxxl«— 150  SACF  SLO 
Rxx2<— 150  SACF  SL1 
Rxx3<— 150  SACF  SL2 
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Figure  D.21.  SL2  Data  Set  Using  Time  Series  Procedure. 


Ti  C  N3 

Figure  D.22.  SL3  Data  Set  Using  Time  Series  Procedure. 

Rxx4«— 150  SACF  SL3 

\ 

Using  the  X-Y  plotting  procedure  in  STATGRAPHICS  the  estimated  auto¬ 
correlation  is  plotted  as  shown  in  Figures  D.23  to  D.26. 

Another  method  for  producing  the  autocorrelation  values  and  plot  them  is  to 
use  the  autocorrelation  procedure  in  the  TIME  SERIES  ANALYSIS  section  from 
the  STATGRAPHICS  Menu.  The  results  are  shown  in  Figures  D.27  to  D.30. 
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Figure  D.27.  Estimated  Autocorrelation  Using  Time  Series  Procedure. 
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Figure  D.28.  Estimated  Autocorrelation  Using  Time  Series  Procedure. 
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Figure  D.29.  Estimated  Autocorrelation  Using  Time  Series  Procedure 
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Figure  D.30.  Estimated  Autocorrelation  Using  Time  Series  Procedure 
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EC  3410  COMPUTER  ASSIGNMENT  [5] 

Compute  power  spectrum  estimates  for  the  data  sets  SLO,  SL1,  SL2,  and  SL3 
using  a  rectangular  window.  In  all  cases  zero-pad  your  data  appropriately  before 
giving  it  to  the  FFT  routine  so  that  you  will  have  computed  enough  points  of  the 
spectrum  to  give  a  smooth  plot.  (Plotting  say  256  points  of  the  spectrum  should 
be  sufficient).  Assume  the  data  has  been  sampled  at  a  spacing  of  T  =  0.1  ms. 
Label  the  frequency  axis  in  hertz.  Since  the  signals  axe  real,  it  is  only  necessary  to 
plot  the  positive  frequencies. 

a.  Use  only  one  16-point  segment  (the  first  16  points)  of  each  data  set. 

b.  Use  32  16-point  segments. 

c.  Use  16  32-point  segments. 

d.  Use  4  128-point  segments. 

What  do  you  make  of  the  results? 


SOLUTION 

1.  The  first  step  is  to  extract  the  first  16  samples  of  each  data  set  and  compute 
the  power  spectrum  estimate  using  the  function  PSE.  In  order  to  obtain  a 
smooth  plot  the  data  sets  are  zero  padded  with  240  points.  (See  Figure  D.31.) 
The  power  spectrum  estimate  will  have  250  points  as  shown  on  the  plots  of 
Figures  D.32  to  D.35. 

2.  By  using  the  function  PSEB  (po./er  spectrum  estimation  with  Bartlett  win¬ 
dow),  we  estimate  the  periodogram  of  the  SLO  data  set  using  30  segments  of 
17  points,  16  segments  of  32  points,  and  4  segments  of  128  points  as  seen  in 
Figure  D.36.  Figure  D.37  is  the  periodogram  plots  for  each  case.  The  above 
process  is  repeated  with  the  data  sets  SL1,  SL2,  and  SL3.  See  Figures  D.38 
to  D.43. 


3.  From  these  results  we  conclude  that  the  first  data  set  has  a  spectrum  like 
white  noise.  The  second  and  third  spectra  are  lowpass  with  bandwidths  of 
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Figure  D.31.  Calculating  the  Power  Spectrum  Estimates  Using  the 
Function  PSE. 


approximately  1  kHz  and  1.2  kHz.  The  fourth  spectrum  is  highpass  with 
bandwidth  of  approximately  1  Khz  Hz.  It  seems  that  the  best  estimate  is 
obtained  by  using  16  segments  of  32  points  each. 
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FREQUENCY  C  KHz 3 


Figure  D.34.  The  PSE  of  SL2  Data  Set  Using  16  Points, 
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Figure  D.35.  The  PSE  of  SL3  Data  Set  Using  16  Points 
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Figure  D.36.  Calculating  Periodograms  of  SLO  Using  Bartlett  Win 
dows. 
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Figure  D.38.  Calculating  Periodograms  of  SLl  Using  Bartlett  Win¬ 
dows. 
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Figure  D.39.  Periodograms  of  Different  Windows 
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Figure  D.40.  Calculating  Periodograms  of  SL2  Using  Bartlett  Win- 
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Figure  D.42.  Calculating  Periodograms  of  SL3  Using  Bartlett  V^in 


EC  4440  COMPUTER  ASSIGNMENT-FILTER  DESIGN  [6] 

This  problem  deals  with  the  design  of  a  lowpass  filter  by  the  windowing 
method.  Use  APL  to  do  this  problem.  Use  STATGRAPHICS  to  generate  3-D 
and  contour  plots. 

1.  Start  with  the  ideal  frequency  response  defined  by  the  following  equation: 


/(u».u,)  =  (1  0-4*)2 

\  0  ,  otherwise  —  n  <  uq, 


U>2  ^  7T 


Construct  an  APL  array  with  its  values.  Note  that  for  compatibility  with  the 
FFT,  the  range  on  uq  and  uq  should  be  zero  to  2  rather  than  — tt  to  +7r.  Use 
a  32  by  32  point  array. 

2.  Find  the  ideal  impulse  response  i(nl,n2)  and  plot  it. 

3.  Choose  one  of  the  1-D  windows  and  plot  it.  Define  a  2-D  rectangular  window 
from  it.  Multiply  t(nl,n2)  by  the  window  to  get  h(nl,n2).  The  window 
should  by  1 1  by  1 1  points  in  size. 

4.  Use  the  Fourier  Transform  (FFT)  of  h(nl,n2 )  to  get  the  filter  frequency  re¬ 
sponse.  Plot  this  using  both  3-D  plotting  and  contour  plotting. 


SOLUTION 

1.  By  using  the  function  ILPFILT  (Ideal  LP  filter)  samples  of  the  freauencv 
response  of  an  ideal  lowpass  filter  is  generated  (32  by  32  point  array).  The 
plots  achieved  by  using  the  surface  plotting  procedure  in  STATGRAPHICS 
can  be  seen  in  Figure  D.44. 

2.  By  using  the  IFFT2D  function  (Inverse  2-D  Fast  Fourier  Transform)  the  ideal 
impulse  response  is  founded  and  plotted  in  Figure  D.45. 

3.  By  using  the  function  HAMMINGR  we  generate  a  2-D  HAMMING  window 
with  rectangular  base.  See  Figure  DAG. 

4.  Figure  D.47  shows  the  product  of  the  HAMMING  window  and  the  ideal  im¬ 
pulse  response. 
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Figure  D.45. 
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Figure  D.46.  Generating  and  Plotting  the  Hamming  Window. 
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Figure  D.47.  The  Product  of  the  Hamming  Window  and  the  Ideal  Im¬ 
pulse  Response. 
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Figure  D.48.  The  Design  Frequency  Response. 
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Figure  D.49.  The  Shifted  Design  Frequency  Response 
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